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ABSTRACT 

Energy released by a small fraction of the baryons in the universe, which condensed 
out while the IGM was cold, dark, and neutral, reheated and reionized it by redshift 
6, exposing other baryons already condensed into dwarf-galaxy minihalos to the glare 
of ionizing radiation. We present the first gas dynamical simulations of the photoe- 
vaporation of cosmological minihalos overtaken by the ionization fronts which swept 
through the IGM during the reionization epoch in the currently-favored ACDM uni- 
verse, including the effects of radiative transfer. These simulations demonstrate the 
phenomenon of I-front trapping inside minihalos, in which the weak, R-type fronts 
which traveled supersonically across the IGM decelerated when they encountered the 
dense, neutral gas inside minihalos, and were thereby transformed into D-type I-fronts, 
preceded by shock waves. For a minihalo with virial temperature below lO^K, the I- 
front gradually burned its way through the minihalo which trapped it, removing all of 
its baryonic gas by causing a supersonic, evaporative wind to blow backwards into the 
IGM, away from the exposed layers of minihalo gas just behind the advancing I-front. 
We describe this process in detail, along with some of its observable consequences, for 
the illustrative case of a minihalo of total mass 10^ Mq, exposed to a distant source of 
ionizing radiation with either a stellar or quasar-like spectrum, after it was overtaken 
at redshift z = 9 by the weak, R-type I-front which ionized the IGM surrounding the 
source. For a source at z = 9 which emits 10^^ ionizing photons per second at 1 Mpc 
(or, equivalently, 10^^ ionizing photons per second at 10 kpc), the photoevaporation 
of this minihalo takes about 100-150 Myrs, depending on the source spectrum, ending 
at about z =7-7.5. 

Such hitherto neglected feedback effects were widespread during the reionization 
epoch. N-body simulations and analytical estimates of halo formation in the ACDM 
model suggest that sub-kpc minihalos such as these, with virial temperatures below 
lO^K, were so common as to cover the sky around larger-mass source halos and possibly 
dominate the absorption of ionizing photons during reionization. This means that 
previous estimates of the number of ionizing photons per H atom required to complete 
reionization which neglected this effect may be too low. Regardless of their effect on 
the progress of reionization, however, the minihalos were so abundant that random 
lines of sight thru the high-z universe should encounter many of them, which suggests 
that it may be possible to observe the processes described here in the absorption 
spectra of distant sources. 

Key words: hydrodynamics — radiative transfer — galaxies: halos — galaxies: high- 
redshift — intergalactic medium — cosmology: theory 
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1 INTRODUCTION 

1.1 The reionization epoch 

Observations of the intergalactic medium (IGM) indicate 
that the universe was reionized by z « 6 but that reion- 
ization began much before this. The absence of a Gunn- 
Peterson (GP) trough in the spectra of high-redshift quasars 
at observed wavelengths Aa(l+2) due to Lya resonance scat- 
tering by H atoms in the IGM at 2 < 6 indicates that the 
IGM was very highly ionized at all 2 < 6 ( Fan ct al. 2000). 
Since the universe recombined at 2 ~ 1000, something must 
have occurred to reionize it by 2 ~ 6. 

The detection of GP troughs in the spectra of SDSS 
quasars at 2 ^ 6 places a lower limit on the mean H I 
density in the IGM at z — 6 large enough to suggest 
that reionization might only just have ended at 2 » 6 
jBecker et al.ll200ll: iDiorgovski et alJl200ll: iFan et al.ll2003l) . 
Limits on the GP effect at 2 < 6, that is, imply a hydro- 
gen neutral fraction significantly smaller than the value at 
z — 6, which can be explained by a photoionized IGM only 
if the ionizing background radiation increased dr amatically 
in a short interval of cosmic time at this epoch jPan et alJ 
L2OO2; Ccn & McDonald 2002). Since reionization is believed 
to have proceeded inhomogeneously, with isolated patches 
of ionized gas which eventually grew to fill all of space, the 
epoch of overlap would have been accompanied by the re- 
quired rapid increase in the ionizing radiation background, 
as the absorption mean free path for this radiation suddenly 
increased. For this reason, the GP detections ai z — 6 and 
limits at 2 < 6 are thought to indicate that the epoch of 
overlap - the end of reionization - occurred at a redshift 
close to 2 = 6. 

The recent discovery by the WMAP satellite of polar- 
ization of the CMB which fluctuates on an angular scale of 
10°, on the other hand, amounts to a detection of a fore- 
ground electron sca ttering optical depth through the IGM 
of res ~ 0.17 ±0.04 jKoeut et al."2003l: ISpergel et a l."2003). 
This suggests that the IGM was either mostly ionized by a 
redshift of 2 ~ 17 ± 5, or else was at least pa rtially ionized 
starting at even higher 2 iSoereel et al.ll2003l) . The GP de- 
tection at 2 = 6 described above argues against the simplest 
interpretation of the WMAP polarization result in which the 
universe was fully reionized by redshift 2 > 15 and stayed 
ionized thereafter. Another argument against that simple 
interpretation is based upon the fact that, after the IGM 
was reheated by its reionization, it would gradually have 
cooled by adiabatic expansion even though it continued to 
be photoionized by background radiation. If reionization fin- 
ished as early as implied by the simplest interpretation of 
the WMAP results, this would have lowered the tempera- 
ture of the IGM at re dshifts 2 4 below the level deduced 
from the Ly a forest (iHui fc Hai man 2003). 

In short, reionization is currently believed to have begun 
by redshift 2 > 15 to accomodate the WMAP results but to 
extend until 2 = 6 to accomodate observations of the GP 
effect and of the Ly a forest at 2 ^ 6. When and how this 
reionization occurred and what its consequences were for the 



subsequent epochs is one of the major unsolved problems of 
cosmology. 

1.2 Ionization fronts in the IGM 

The neutral, opaque IGM out of which the first bound ob- 
jects condensed was dramatically reheated and reionized at 
some time between a redshift 2 ~ 30 and 2 ~ 6 by the 
radiation released by some of these objects. When the first 
sources turned on, they ionized their surroundings by propa- 
gating weak, R-type ionization fronts which moved outward 
supersonically with respect to both the neutral gas ahead 
of and the ionized gas behind the front, racing ahead of 
the hydrodynamical respons e of the IG M, as first described 
by^Shaoiro (1986) and Shap iro fc Girou x (1987). These au- 
thors solved the problem of the time-varying radius of a 
spherical I-front which surrounds an isolated source in a 
cosmologically-expanding IGM analytically, generalizing the 
I-front jump condition and radiative transfer to cosmolog- 
ical conditions. They applied these solutions to determine 
when and how fast the I-fronts surrounding isolated sources 
would grow to overlap and, thereby, complete the reioniza- 
tion of the universe^. The effect of density inhomogeneity 
on the rate of the I-front propagation was described by a 
mean "clumping factor" C = (n^)/{n)^ > 1, which slowed 
the I-fronts by increasing the average recombination rate 
per atom inside clumps. This suffices to describe the rate 
of I-front propagation as long as the clumps are either not 
self-shielding or, if so, only absorb a small fraction of the 
ionizing photons emitted by the central source. 

Numerical radiative transfer methods are currently un- 
der development to solve this problem in 3-D for the inho- 
mogeneous density distribution which arises as cosmic struc- 
ture forms, so far primarily limited to the transfer of ioniz- 
ing radiation through pre-computed density fields generated 
by cosmological simulations which do not include radiative 
transfer, so there is no back-reaction of the radiative transfer 
calculati on on the density field or any aspect of the gas dy- 
namics azoumov fc Sc ott '1999*: lAbel. Norman fc Madaul 
1999l:ICiardi et all200lHNakamoto Umemura fc Susal200l|: 
Sokasian et alJ I2001I: ICenI 120021: iHaves fc NormanI l2002l : 
Ciardi et al.ll2003ll . Approximate approaches which incor- 
porate some of the important effects of radiative trans- 
port during reionization within the context of cos mological 
gas d ynamics simulation have also been developed (jCnedi nl 
200dl: lGnedin fc Abel20(5ll: iRicotti. Gnedin fc Shull2002a^ . 
These recent attempts to model inhomogeneous reioniza- 
tion numerically are generally handicapped by their lim- 
ited spatial resolution, which prevents them from fully re- 
solving the most important (sub-kpc-sized) density inhomo- 
geneities. The questions of what dynamical effect the I-front 

^ This provided the first indication that quasars alone did not 
reionize the universe; the observed quasar luminosity function 
extrapolated to higher redshift was not sufficient to do so in time 
to satisfy the GP constraint, suggesting that some other source 
was responsible, such as starhght. 
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had on the density inhomogeneities it encountered and how 
the presence of these inhomogeneities affected the f-fronts 
and reionization, therefore, require further analysis. 

Towards this end, we have developed a radiation- 
hydrodynamics code which incorporates radiative transfer 
and focused our attention on properly resolving this small- 
scale structure. Here we shall present in some detail our re- 
sults of the first radiation-hydrodynamical simulations of the 
back-reaction of a cosmological I-front on a gravitationally- 
bound density inhomogeneity it encounters - a lO'^M© dwarf 
galaxy minihalo - during reionization, along with some gen- 
eral considerations which apply to other halos as well. A 
second paper to follow will present the results for a much 
wider range of cases, with different halo masses and I-front 
encounter epochs, to quantify just how the process described 
here varies with these properties. As it turns out, the photo- 
evaporation of sub-kpc-sized objects like these may be the 
dominant process by which ionizing photons are absorbed 
during reionization, so this problem is of critical impor- 
tance. In addition, observations of the absorption spectra of 
high redshift sources like those which reionized the universe 
should reveal the presence of these photoevaporative flows 
and provide a useful diagnostic of the reionization process. 

1.3 The photoevaporation of dwarf galaxy 
minihaloes overtaken by cosmological 
ionization fronts 

The effect which small-scale dumpiness had on reioniza- 
tion depended upon the sizes, densities, and spatial dis- 
tribution of the clumps overtaken by the I-fronts dur- 
ing reionization. For the currently-favored ACDM model 
(fio = 1 - Ao = 0.3, h = 0.7, D.th^ = 0.02, primor- 
dial power spectrum index Up — 1; COBE-normalized) , 
which is the model we adopt throughout this paper, the 
universe at 2 > 6 was already filled with dwarf galaxies 
capable of trapping a piece of the global, intergalactic I- 
fronts which reionized the universe and photoevaporating 
their gaseous baryons back into the IGM. Prior to their en- 
counter with these I-fronts, "minihalos" with Tvir ^ lO^K 
were neutral and optically thick to hydrogen ionizing radi- 
ation, as long as their total mass exceeded the Jeans mass 
Mj in the unperturbed background IGM prior to reioniza- 
tion (see §|lll, as was required to enable baryons to collapse 
into the halo along with dark matter. Their "Stromgren 
numbers" Ls = 2i?haio/^s, the ratio of a halo's diameter 
to its Stromgren length is inside the halo (the length of 
a column of gas within which the unshielded arrival rate 
of ionizing photons just balances the total recombination 
rate), were large. For a uniform gas of H density nn, lo- 
cated a distance rMpc (in Mpc) from a UV source emitting 
A^ph,56 ionizing photons (in units of 10^® s~^), the Stromgren 
length is only £3 « (100pc)(iVph,56/?-Mpc)(wH/0.1 cm-^)-^ 
so Ls ^ 1 for most of the range of halo masses and sources 
of interest, as we will discuss in fi \'2.3\ In that case, the inter- 
galactic, weak, R-type I-front which entered each minihalo 
during reionization would have decelerated to about twice 



the sound speed of the ionized gas before it could exit the 
other side, thereby transforming itself into a D-type front, 
preceded by a shock. Typically, the side facing the source 
would then have expelled a supersonic wind backwards to- 
ward the source, which shocked the surrounding IGM as the 
minihalo photoevaporated. 

The importance of this photoevaporation process 
has long been recognized in the study of interstellar 
clouds exposed to ionizing starlig ht (|Oortj£.&oitzer 
19551: ISpitzeil ITqtI iBertoldil llOSa iBertoldi fc McKee 
199C : jLefioch & Lazareflf 'l994'; 'Bertoldi & Drains 
I995I: iLizano ot al.^ 1996, ; Gorti fc HoUenbach 200 JT 
Radiation-hydrodynamical simulations were performed 
in 2-D of a stellar I-front overtaking a clump inside 
a molecular cloud jKlein. Sandford fc Whitaked Il982l: 
ISan dford. Whitaker fc KleinI \l98^ . More recentlv. 2-D 
simulations for the case of a circumstellar clo ud ionized by a 
single nearby star have also been performed iMellema et alJ 
Il998h . In the cosmological context, however, the importance 
of this process has only recently been fully appreciated. 

In proposing the expanding minihalo model to ex- 

plain Lyman a for est ["LP") quasar absorption lines, 

iBond. Szalav fc Silk! lll988r discussed how gas originally 
confined by the gravity of dark minihalos in the CDM model 
with virial temperatures below 10^ K would have been ex- 
pelled by pressure forces if photoionization by ionizing back- 
ground radiation suddenly heated all the gas to an isother- 
mal condition at T « 10* K, a correct description only in the 
optically-thin limit. The first realistic discussion of the pho- 
toevaporation of such minihaloes by cosmological I-fronts, 
including the first r adiation-hydrodynamical s i mulations o f 
this process, was bv lShapiro. Raea fc Mellemal Jl997l. Il998ll . 
This work demonstrated the importance of taking proper 
account of optical depth and self-shielding simultaneously 
with hydrodynamics, in order to derive the correct dynami- 
cal fiow, ionizat ion structure and t he res ulting observational 
characteristics. iBarkana fc Loebl lll999ll subsequently con- 
firmed the relative importance of this process as a feedback 
effect on dwarf galaxy minihaloes during reionization, using 
static models of uniformly-illuminated spherical clouds in 
thermal and ionization equilibrium, without accounting for 
gas dynamics. They took H atom self-shielding into account 
and assumed that all gas which is instantaneously heated 
above the minihalo virial temperature must be evaporated. 
They concluded that 50%-90% of the gas in gravitationally- 
bound objects when reionization occurred should have been 
evaporated. However, by neglecting the time-dependent gas 
dynamical nature of this phenomenon, their model failed to 
capture some essential physics, leading to incorrect results, 
as we will show below. 

Further results of radiation-hydrodynamical simula- 
tions of a cosmological minihalo overrun by a weak, R- 
type I-front in t he surrounding IG M were summa rized in 
IShapiro fc Rag3 (ipOOa b, 200% and lShapird J200ll) . These 
results were preliminary versions of the final, more advanced 
results which we shall present in this paper in some detail. 
Before we do, it is worth summarizing why we believe that 
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the particular problem simulated here is not only generic to 
the minihalos which formed prior to the end of the reioniza- 
tion epoch, but likely to have affected the global outcome of 
cosmic reionization, as well. 

According to the currently prevalent paradigm of cos- 
mic structure formation - the flat, cold-dark matter universe 
with cosmological constant (ACDM) - cosmological struc- 
tures formed hierarchically, with small mass haloes form- 
ing first and merging over time to form ever larger halos. 
The first baryonic structures to emerge in this way were 
the "minihalos" which began forming at z > 20 — 30, with 
masses in the range from about IO'^A/q to 10* Mq, where the 
minimum is set by the Jeans mass in the cold and neutral 
IGM prior to its reionization. The term "minihalo" used here 
specifically refers to halos at any epoch with virial temper- 
atures Tvir ^ 10* K. In such halos, the thermal coUisional 
ionization rate alone was not sufficient to ionize the H and 
He atoms significantly, so the minihalo gas was almost en- 
tirely neutral, in the absence of ionizing radiation. Moreover, 
coUisional line excitation of H atoms is exponentially sup- 
pressed below 10* K and, thus, a purely atomic gas of H and 
He could not radiatively cool effectively below that temper- 
ature. Such minihalo gas was able to cool below 10* K only 
if H2 molecules formed in sufficient abundance to cool it 
by rotational- vibrational line excitation, to temperatures as 
low as T < 100 K. This H2 formation, catalyzed by the tiny 
residual ionized fraction, was essential if the minihalo gas 
was to cool and, thereby, compress relative to the dark mat- 
ter, enough to become self-gravitating and form stars. Sim- 
ulations suggest that minihalos of mass M>lifMQ did just 
that, thereby forming the first me tal- free ("Pop. HI") stars 
of mass Mt > IQOMp) at z < 30 l|Bromm. Coppi fc LarsonI 

), perhaps leading to the 



12002: lAbel. Brvan fc Norman 
first miniquasars, as well. 

While this process of star formation in minihalos is im- 
portant as the origin of the first sources of ionizing radiation 
in the universe, minihalos are not likely to have been the 
dominant source of cosmic reionization. The H2 formation 
required for minihalos to form stars would have been sup- 
pressed easily by photodissociation in the Lyman- Werner 
bands by the background of UV radiation created by the 
very first minihalos which formed stars, when the ioniz- 
ing radiatio n background was still much too low to cause 
reionization ijHaiman. Abel fc ReejEoOff) . If so, most mini- 
halos must not have formed stars and must, instead, have 
remained "sterile" , dark-matter-dominated halos filled with 
stable, neutral atomic gas at the halo virial temperature 
Tvir < 10* K. In that case, when some other source of ra- 
diation ionized the universe and drove an I-front into the 
minihalo, which photoheated the gas to T > 10* K, mini- 
halo gravity would have been unable to confine the ion- 
ized gas and, eventually, all of it was expelled. A caveat 
to this description above of the negative feedback on the 
H2 in minihalos due to the first UV sources is the possi- 
bility that the presence of a weak X-ray background might 
have counterbalanced this H2 photodissociation by enhanc- 
ing the H2 formation rate in places where X-ray ioniza- 



tion increased the sma ll ionized fraction inside the halos 
jHaiman. Abel fc Reed |2000).^ Apart from such caveats, 
however, the primary source of reionizing photons is believed 
to have been the halos with Tvir ^ 10* K. For such halos, 
radiative cooling by H atom line excitation alone was effi- 
cient enough to compress the gas to enable self-gravity and 
self-shielding, leading to star formation. For the purposes 
of this paper, we shall henceforth accept the premise that 
the minihalos were predominantly "sterile" targets of source 
halos with Tvir > 10* K, of higher mass. 

In that case, it has been shown, source halos would 
generally have found their skies covered by the miniha- 
los between them and neighboring source halos, until the 
ionization fronts created by those sources overtook the 
minihalos and photoevaporated them, possibly consuming 
a significant fraction of the ionizing photons released by 
the sources in the process (plairnan, Abel fc Madau 20M 
[Shapird l200ll . iBarkana fc Loebil2002il and Ishapiro et'alJ 
l2003h . According to these authors, the phenomenon of 
minihalo photoevaporation during reion ization identified by 
IShapiro. Raea fc Mellernal lll997l Il998i) not only was quite 
common and had an important negative feedback effect on 
the minihalos, it also had an important effect on reioniza- 
tion, itself, by screening the ionizing sources and increas- 
ing the number of ionizing photons required per baryon to 
complete reionization. While these results demonstrate the 
importance of minihalos to the cosmic reionization story, we 
will show that a correct quantitative estimate of these ef- 
fects requires detailed gas dynamics and radiative transfer 
simulations. In addition, as pointed out by Shapiro (200^, 
there is a very large number of minihalos per unit redshift 
along an arbitrary line-of-sight (LOS) through the universe 
at z ^ 6, so the process of their photoevaporation is po- 
tentially observable in the absorption spectra of very high-z 
sources. 

This paper is organized as follows. In § |21 we discuss 
our minihalo model and the statistics of minihalos in the 
ACDM universe, as well as the analytical expectations for 
how the photoevaporation of cosmological minihalos pro- 
ceeds and what its outcome is. In § |3 we outline the basic 
equations used to model the gas dynamics of photoevapora- 
tion, the microphysical processes considered and the initial 
conditions for our simulations. In §2] we discuss our numeri- 
cal method and tests, and the parameters of our simulations. 
In §1^1 we describe the photoevaporation process in detail for 
the illustrative case of a minihalo of 10^ Mq at 2 = 9 exposed 



^ The net effect of these early X-rays on the H2 in minihalos 
may have beem more negative than positive, however, if they 
raised the entropy floor of the IGM from which the minihalos 
condensed JOh &: Haimanl2003l) . Another study suggests that the 
near environment of early-star forming halos might also have ex- 
perienced some positive feedback on the H2, which partially off- 
set the negative feedback due to H2 dissociation by UV sources, 
although not necessarily insid e the pre-existing minihalos (e.g. 
iRicotti. Gnedin fc Shull2002bD : this study did not include an X- 
ray background. 
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to ionizing radiation from one of three types of sources: (1) 
starlight with a 50, 000 K black-body spectrum (hereafter, 
referred to as "BB 5e4"), representative of massive Pop. II 
stars, (2) starlight with a 10^ K black-body spectrum (here- 
after, "BB le5"), as expected for massive Pop. Ill stars, 
and (3) QSO-like, with a power-law spectrum F^, oc u~^'^ 
{u > vh ) (hereafter, "QSO"). We show our results for the 
structure of the global I-front propagating in the IGM dur- 
ing its weak, R-type phase (S 15.11 . the encounter between 
this global I-front and a minihalo along its path and the 
penetration of the minihalo by the front in its weak, R-type 
phase (§ 5.2), the trapping of the I-front inside the mini- 
halo and its transition to D-type (§ 5.3), the structure of 
the photoevaporative flow during the D-type phase I5.4II . 
the evolution of temperature structure (§ 5.5) and I-front 
location and speed (§ 5.6), the evaporation times (§ 15.711 . 
the consumption of ionizing photons (S 15. St . and some ob- 
servational diagnostics 15.91 . Finally, in § 1^1 we give our 
summary and conclusions. 



2 ANALYTICAL CONSIDERATIONS 

2.1 Cosmological minihalos during reionization in 
a ACDM Universe 

2.1.1 Minihalo model 

Our initial minihalo model consists of a virialized region in 
hydrostatic equilibrium embedded in an exterior region of 
cosmological infall. For the former, we have previously de- 
veloped an analytical model for the postcoUapse equilibrium 
structure of virialized objects that condense out of a cosmo- 
logical background universe, eit her matter-dominated or flat 
with a cosmological constant JShapiro. Iliev fc Ragalll999t 
llliev fc Shapiro! l200l[) . This Truncated Isothermal Sphere, 
or TIS, model assumes that cosmological halos form from 
the collapse and virialization of "top-hat" density perturba- 
tions and are spherical, isotropic, and isothermal. This leads 
to a unique, nonsingular TIS, a particular solution of the 
Lane-Emden equation (suitably modifled when A 7^ 0). The 
size rt and velocity dispersion ov are unique functions of the 
mass M and formation redshift Zcoii of the object for a given 
background universe. The TIS profile fiattens to a constant 
central value, po, which is roughly proportional to the criti- 
cal density of the universe at the epoch of collapse, Pc(zcoii), 
with a small core radius ro « rf/30 [where ay — A-kGpqTq 
and rp = rKing/3 fo r the "King radius" rKing, defined by 
iBinnev fc Tremaine I lll98i1) . p. 228]. While the TIS model 
does not produce the central cusp at very small radii in the 
dark matter density profile of halos predicted by numerical 
CDM N-body simulations, it does reproduce many of the av- 
erage properties of these halos remarkably well, suggesting 
that it is a useful approximatio n for halos which result from 
more realistic initial condit ions dShapiro. Iliev fc Ragall999l: 
llliev fc Sh apiro'2001'. 2002). In particular, the TIS mass pro- 
file agrees well with the fit by NFW to numerical simulations 
(i.e. fractional deviation <20%) at all radii outside of a few 



TIS core radii (i.e. outside a King radius or so). Our appli- 
cation here of the TIS model is not sensitive to the small 
difference between the TIS and NFW dark-matter density 
profiles at small radii. In both cases, the baryonic profile is 
nonsingular. 

The TIS halo is uniquely specified by the central den- 
sity Po and core radius ro. Since the central density is pro- 
portional to pc(2coii), it is also possible to specify the pro- 
file by the two parameters, total mass and collapse redshift, 
(Af, Zcoii), which is equivalent to specifying the pair (ro, po). 
This makes the TIS model extremely useful in combination 
with other analytical methods which determine the typical 
epochs of collapse for halos of different masses, such as the 
well-known Press-Schechter approximation. We note that it 
is customary elsewhere to define the total mass of profiles 
used to model CDM halos as M200, the mass inside a sphere 
of radius r2oo with a mean density which is 200pc(z) at some 
redshift z. For the sake of comparison with the TIS profile, 
if we were to fix M200 and r2oo for halos of different pro- 
files, which amounts to fixing M and Zcoii for the TIS halo, 
then M = I.I67M200 = 772.6porg, po = 18, 000pc(2coii) and 
r2oo = 24.2ro if 2 = Zcoii- 

As applied here to model minihalos, the TIS is em- 
bedded in a self-similar, spherical, cosmological infall, us- 
ing the solution for the latter i n an Einstein-de Sitter uni- 
verse derived bv iBertschingeJ 1^8^, generalized by us to 
apply it to a low-density background universe at the early 
times considered here (Il iev fc Shapiro 2001, Appendi x A). 
As discussed in detail in IShapiro. Iliev fc Raeal ^ll999^ and 
llliev fc Shapiro! i200 J) , the truncation radius of the TIS so- 
lution coincides almost exactly with the location of the shock 
in the self-similar infall solution, allowing a natural match of 
the two models. The procedure which we used to accomplish 
that is discussed in detail in the Appendix. 

For the calculations in this paper, our fiducial case is 
a minihalo of mass M — 10^ M0 which collapses at z = 
9. According to the TIS model, this minihalo has a virial 
radius rt — 0.76 kpc, a virial temperature Tvir ~ 4000 K, 
and a DM velocity dispersion ay ~ 5.2 km s"'^. In units 
of these fiducial values, we define M7 = M/{W''Mq) and 
(1 + z)io = (1 -f z)/10. We shall express temperature in 
terms of r4 = r/(10''K). 

2.1.2 Statistical properties of minihalos 

This TIS model for the internal structure of an individual 
minihalo of a given mass M which forms at a given collapse 
epoch Zcoii can be combined with the Press-Schechter (PS) 
approximation to determine the statistical expectations for 
the properties of halos of different masses at each epoch 
which result from the Gaussian-random initial fluctuations 
of the ACDM universe. 

In Figure^ we plot the curves in the halo M— Zcoii plane 
which correspond to halos which collapse from v-a fluctua- 
tions for v = 1,2, and 3. Here v = 5crit/cr(M), a{M) is the 
rms fluctuation in the dark matter density fleld, filtered on 
mass scale M, as predicted by linear theory, and Sait is the 
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Figure 1. How common are minihalos at high redshift? Plot of 
minihalo mass vs. redshift. Shaded region represents minihalos: 
halos with virial temperature Tvir < 10* K. Dotted curves repre- 
sent V — a fluctuations with v = 1, 2 and 3, predicted by the PS 
approximation, as labeled. 



overdensity of a top-hat perturbation given by extrapolat- 
ing the hnear solution to the time of infinite collapse in the 
exact nonlinear solution. The value v = 1 corresponds to 
the most common halos, while high values of correspond 
to rare density peaks. 

The minihalos span a mass range from Afmin to Mmax 
which varies with redshift, with Afmin close to the Jeans 
mass of the uncoUapsed IGM prior to reionization, 



Mj 



while 



5.7 X IG'M, 







0.15 



-1/2 



0.02 



-3/5 



x(l+^)?/^ 



3.95 X 10^ Mq 



0.15 



-1/2 



(1 + -) 



-3/2 



(1) 



(2) 



is the mass for which Tvir = 10* K according to the TIS 
model (Iliev & Shapiro 2001). This maximum minihalo mass 
at each epoch (i.e. the curve of constant Tvir = 10* K) 
is shown in Figure Figure demonstrates that, below 
redshift z ~ 10, all minihalos result from fluctuations with 
I' < 2, so they are quite typical and fairly common. In- 
tegrating over the entire mass function of minihalos from 
Mmin to Mmax, wc fiud that the total collapsed mass frac- 
tion in minihalos is (36%,28%,10%,3%) at z = (6,9, 15,20), 
respectively. 



2.2 Ionizing flux 

We shall assume that each minihalo is photoevaporated by 
an isotropic point source located at a proper distance r from 
the minihalo centre, which emits A'^ph H ionizing photons per 
second with luminosity L^(erg s~^Hz~^), so that 



iVph = 



hv 



(3) 



The unattenuated flux at the location of the minihalo cen- 
tre is just = Liy/(47rr^), since r <^ c/H{z), the hori- 
zon size during minihalo photoevaporation. We shall, hence- 
forth, parametrize the flux to which a given halo is exposed 
in terms of the dimensionless frequency-integrated photon 
flux, Fo = A'ph.se/fMpc- For example, a source which emits 
10^^ (or lO^'^) ionizing photons per second at a distance of 
1 Mpc (or 10 kpc) corresponds to Fq = 1, or, in physical 
units, a flux of 8.356 x lO^cm'^^s"^. What values of _Fo are 
relevant to the encounter between a minihalo and an I-front 
during cosmic reionization? 

To get a sense of how large this fiducial fiux is in terms 
of an equivalent isotropic radiation background, we can esti- 
mate the mean ionizing photon intensity J-,(cm~^s~^ster~^) 
required to irradiate an atom on the minihalo surface with 
the same fiux, by equating this fiux to 2nJj. For a QSO- 
like ionizing background with a power-law energy spectrum, 
Ji, cx a > 1, for example, the integrated photon 

intensity is J-y — 1.509 x 10^J_2i(l + a)~^, where J_2i 
is the mean intensity at the H Lyman limit in units of 



10"^^ergcm"^s"^Hz" 



^ster ^ . The equivalent dimensionless 



fiux parameter is then given by 
fo.oqu = 1.13(1 -Hq)-'J_21, 



(4) 



so Fo = 1 is roughly equivalent to J_2i ^ 1. 

We can estimate the appropriate range to consider 
for Fo as follows. Suppose the luminosity of a dwarf 
galaxy is prod u ced b y a super-massive star cluster (SSC). 
iTan fc McKed i200 estimate the ionizing luminosity of 
such an SSC, using the observed stellar mass function in the 
R136 cluster in 30 Doradus as input to the STARBURST99 



code of lLeitherer et all il9& 
100): 

A'ph,52 = 3.94M. 



finding (for 0.1 < A'/,/Afg 



(5) 



where Af« is the total mass in stars, Mp = M / (10^ Mq) and 
iVph,q = iVph/(10'' s ~^) is the fiux o f ionizing photons. 

Alternatively, iMadau fc ReesI (l2000l) have estimated 
that, for a standard Salpeter IMF, approximately 3000 — 
4000 ionizing photons are emitted for every stellar baryon. 
Hence, if the galaxy gas mass turned into stars is Af», the 
total number of ionizing photons emitted during the lifetime 
of those stars is 



Assuming an average lifetime of rufctin 

A'ph,52 ~ 1.6Af,,6, 



(6) 

10^ yrs, we obtain 
(7) 
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similar to the result in equation (j^. On the other hand, 
iBromm. Kudritzki fc Loebl i200 J) estimate that the number 



of photons produced per baryon by Pop. Ill stars can be 
about 10-20 times higher than this. 

If we let /i, — Qi,/Qo ~ 0.136 be the cosmic mean bary- 
onic mass fraction and /* be the fraction of these baryons 
which are converte d into stars in t h e dwa rf galaxy halo, then 
M. = /./bMhaio- iTan fc McKe3 ||200J) find /. = 0.5 for 
their SSC. If the target minihalo mass Mtargct ^ Mhaio, 
the mean distance between source and halo is r ~ dscp/2, 
while if Mtargct ~ Mhaio, it is r = dscp- For source halos 



of mass 10 Mq at z = 9, d. 



i 50 kpc (Shaoiro 2001), so 
ryipc ~ (0.5 — 1) X dsep.Mpc ~ 0.025 — 0.05. The character- 
istic dimensionless flux Fq for such sources at this epoch is, 
therefore, given by 



Fo 



= (1 - 4) X 3.94 X 10-^ ^ 1 - 



' Mpc 



(P 

"sop, Mpc 



-4,(8) 



if we adopt the lum inosity per gas mass estimated by 
iTan fc McKe3 i200 J) . For halos in the dwarf galaxy mass 
range, the PS approximation for the halo mass function 
in ACDM yields drihaio/dM oc l|shapiro 2001), so 

the mean separation of dwarf galaxy sources of mass M is 
dscp = (Mdrihaio/dM)"^/^ oc A/^/^. If the mass-to-light ra- 
tio Ad/L is independent of halo mass, then the flux incident 
on a minihalo due to one of these halos scales only weakly 
with source halo mass, _Fo oc M^^^, up to the galaxy mass 
scale above which the PS mass function cuts off exponen- 
tially, at which point Fq drops more rapidly with increasing 
mass. 

Finally, we can get a simple, direct estimate of _Fo based 
on the mean flux required to reionize the universe, as follows. 
Suppose the universe was filled with sources that supplied ^ 
ionizing photons per H atom steadily over the time between 
the source turn-on epoch Zon and the epoch of reionization 
overlap at Zov ■ Assuming that all photons remained energetic 
enough to ionize H atoms despite their continuous redshift- 
ing, and neglecting the absorption of these photons during 
the extended reionization epoch, the mean comoving num- 
ber density of ionizing photons at any redshift in the range 
Zon ^ z ^ Zov would be just 

t(z) - ton 



■to 



(9) 



where t is the age of the universe at each redshift, as labelled, 
and uh here refers to the mean IGM H density. The mean 
photon intensity of this isotropic background of ionizing 
radiation would then be given by 



_ c _ _ gnnc t{z) - ton 

47r 47r tov — ton 



(10) 



This corresponds to an effective photon flux F^s given by 



^UhC t{z) - ton 



2 ^ov 

or a dimensionless flux 
-Fo.cff =2,i(l + z) 



0.02 



t{z) 



(11) 



(12) 



For ton <S tov the quantity in square brackets can be replaced 
by t{z)/tov to yield an upper limit. 



-Fo.cff <3^{l + z[ 



10 ^ 002" 



1+Z 



\l+ZovJ 



or, if Zov = 6 and Q,bh^ = 0.02, 



i^o,off< 1.8^(1 -Hz 



,3/2 



(13) 



(14) 



This suggests that _Fo '~ 1 is a reasonable fiducial choice for 
the flux to which minihalos were exposed during reionization 
before overlap. In this paper we present results for Fo — 1 
only. We defer the study of the effects of varying _Fo to the 
companion paper. 



2.3 Minihalo Stromgren numbers and I-front 
trapping 

2.3.1 I-front trapping 

For a spherical H II region in a uniform interstellar gas, the 
early expansion phase of the weak, R-type I-front ends when 
the I-front decelerates to twice the sound speed of the ionized 
gas, at which point it makes the transi tion from R-c ritical 
to D-critical, preceded by a shock front llSpitzerl|l97d) . This 
transition occurs when the radius of the I-front approaches 
that of the Stromgren sphere. By analogy, we expect our in- 
tergalactic, weak, R-type I-front to make a similar transition 
to D-type after it enters a minihalo when it slows to the R- 
critical speed as it travels up the density gradient inside the 
minihalo, a distance comparable to the Stromgren length 
along its path through the halo. In order to "trap" the I- 
front, therefore, the minihalo size must exceed this length. 
We can estimate the "trapping" condition as follows. 

The Stromgren length £s{r) at impact parameter r is 
given by 



F = 



(15) 



i.e. by balancing the number of recombinations with the 
number of ionizing photons arriving along a given LOS.'^ 
The integration in equation 1151 is done along the LOS, nn 
is the H number density, rie is the electron density, F is the 
flux of ionizing photons, is the Case B recombination 
coefficient for hydrogen, and we have assumed for simplicity 
that only the H in the gas is ionized. We then define the 
"Stromgren number" for the halo as Ls = 2rt/^s(0), where 
rt is the radius of the halo and ^s(O) is the Stromgren length 
for zero impact parameter. If Ls > 1, then the halo is able to 
trap the I-front, while if Ls < 1, the halo would be ionized 
quickly by the passage of the weak, R-type I-front across 

We have adapted cylindrical coordinates (r, x) in which the 
X-axis is the line between the source and the minihalo center, 
the axis of symmetry of this problem. In writing equation 1151 . 
we have assumed that the distance between the source and the 
minihalo is much greater than the size of the minihalo, so we can 
approximate the LOS as parallel to the x-axis. 
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Figure 2. Inverse Stromgren Layer. The Stromgren length along 
the axis, /s(0), in units of the halo diameter 2rt and ISL mass 
^ISL in units of the total mass Mtot vs. Xs/XSO: where xs = 
^HPq^o)' XS,0 is XS for our fiducial halo and Fq = 1. Sym- 
bols indicate the (roughly) corresponding simulation results for 
Misjj/Mtot (star) and £s(0)/(2rf) (square) at the moment of I- 
front transition from R-type to D-type (see § 15. 2i . Dotted line 
indicates the value xs.crit for which Lg = 1. We assume = 2. 



it, which will not slow down enough to be trapped. In the 
latter case, if Ls <^ 1, the halo gas would be uniformly 
photoheated to T 10* K before any mass motions could 
occur, and the pressure gradient would blow the gas apart 
isotropically, long after the 1-front had exited the halo. 
The threshold condition Ls = 1 is equivalent to 



F : 



jPoro 



(16) 



where p = p/po, C = '"/''o, ft ~ Slb/rio is the baryonic mass 
fraction and fiHinn is the mean gas mass per H atom. The 
dependences on the source fiux and on the minihalo mass 
and size at different redshift are combined if we define the 
parameter 



XS 



= Mtot/{W^M@) 
~^ . We have used 



Po^o 



where Mr 
10 



30, 



0.15 



and 
the 



Xs,o 
fact 



-5/3 



that 



(17) 



4.1 
I 



2 J^* fP{Qd(^ — 3.49 and the scalings of po and ro with the 
mass of the halo and its red shift of collapse accord ing to the 
TIS model at high redshift Jlliev fc Shapiroll200ll) . There is 
a direct correspondence between Lg and xs, and the condi- 
tion Ls = 1 is equivalent to xs ~ Xs.crit where 



Xs.crit = 3.3 X 10 T4 ' ( 



-2 3 -1 
g cm s , 



0.15 



(18) 



(where T4 refers to the temperature of the ionized gas) while 
Ls > I and Ls < I correspond to xs > Xs.crit and xs < 
Xs.crit, respectively. We plot £s{0) vs. xs in Figured For our 
fiducial case, xs = XS,o <C Xs.crit, so we expect the I-front 
to be easily trapped by the minihalo in this case. In general, 
the critical minihalo mass which is just large enough to trap 
the I-front by making is = 1 is given by setting xs = Xs.crit 
in equation 1171 to yield 



Mr 



= 1.8 X 10" 



(l + ^)io] 



-15 rpQ/A 

-'4 



0.15 



0.02 



(19) 



For Fo < 10, even the smallest minihalos which formed at 
the last possible moment before the I-front overtook them 
at the end of the reionization epoch at Zov = 6 would have 
been able to trap the I-front. 

For a uniform halo with gas number density {uh), equa- 
tion 1151 1 reduces to 

(20) 

In terms of the I-front trapping condition (Ls = 1), a TIS 
halo with mean density pTis is equivalent to a uniform halo 
of the same size but with mean density given by 

P = /9o(^^^ =0.24po = 34pTis, (21) 

where pTis is the mean density of the TIS. Hence, a uniform 
halo of the same size as the TIS minihalo which is just large 
enough to trap the I-front will do so only if it is 34 times 
more massive than that TIS halo. 

One important implication of this result is that it would 
be very incorrect to use the volume-average mean density 
of a minihalo instead of this much larger density in equa- 
tion H21|l to determine whether a minihalo of a given mass 
will trap the I-front. In particular, if we replace the realis- 
tic TIS profile of the centrally-concentrated minihalo by a 
uniform-density minihalo with the same mass and radius, 
instead, the Stromgren number, Ls^uniform, of this uniform 
halo will be very much smaller than Ls for the TIS minihalo, 
according to 



Ls.uniform = 8.9 X 10"* (xS.crit /xs) = 0.72M7/^(1 + z) 

i 2\ 2 /r^ l2\ -1/3 



xFo-^T-^/* 



(22) 



0.02 J \ 0.15 

Hence, when a TIS minihalo of a given mass is just 
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Figure 3. The ISL shape inside minihalo (halo boundary indi- 
cated by the solid line), according to the ISL approximation for 
x/xsfi = 0.01, 0.1, 1, 10, 100, and 250 (dashed-line curves from 
left to right), xq is the position of the centre of the halo along the 
X-axis. Both the impact parameter r and axis coordinate x — xq 
are in units of the halo radius rt- Source is located outside the 
box, to the left, at a:: = 0, where xg ^ rt- We assume T4 = 2. 

marginally able to trap the I-front (i.e. xs = Xs.crft), a uni- 
form halo with the same mass and radius has a Stromgren 
number Ls ~ 8.9 x 10~*, far below the minimum required to 
trap the I-front! This means that the phenomenon of I-front 
trapping can be missed entirely by numerical simulations of 
reionization which do not fully resolve the internal structure 
of individual minihalos, even when their resolution in mass 
and length are, in principle, high enough to identify mini- 
halos of the correct mass and mean density. For this reason, 
the simulations we shall report here are high enough in reso- 
lution to guarantee that the internal structure of individual 
minihalos is fully resolved. 

2.3.2 Inverse Stromgren Layer (ISL) in minihalo 

We define the Inverse Stromgren Layer as the region in- 
side the minihalo which is the HII region predicted by the 
Stromgren approximation. This region is bounded by the 
minihalo surface nearest to the source and by the surface de- 
fined at each impact parameter r by ^s(r) in equation 11511 . 
This ISL provides an estimate of the minihalo region which 
will be ionized first, during the weak, R-type I-front phase, 
and of the approximate location of the surface inside the 
minihalo where the transition to D-type will occur as the R- 
type front approaches it. In Figure |5| we show the fraction 
of the halo mass which is located inside the ISL, Mish/Mtot, 
vs. xs, given by 

MisL = / 2-Krdr / p{R)dX. (23) 
Jo Jo 

Here R = {r^ + [X - (r? - r2)i/2]2|i/2^ FigurcHwe also 
plot = ^s(0)/(2rt) (the Stromgren length for r = in 
units of the halo diameter) vs. xs, and we indicate the cor- 
responding simulation results at the approximate moment 



of transition of the I-front from R-type to D-type, to be dis- 
cussed in i) 15.21 At Xs = Xs.cnt, the ISL includes the entire 
minihalo, thus MisL/Mtot = 1 and IsiQ) = 2rt. In Figure El 
we plot the shape of the ISL surface for Xs/xsfi = 0.01, 0.1, 
1, 10, 100, and 250. We see that, even for low values of the 
parameter xs (equivalent to low external flux levels) a sig- 
niflcant fraction of the halo volume becomes ionized quickly 
due to the relatively low density in the halo outskirts. How- 
ever, from Figure|21we see that this volume still corresponds 
to a small fraction of the total mass (e.g. for Xs/xsfi < 0.01, 
MisL < 10%). As Xs approaches xs.crit, the neutral gas frac- 
tion shielded by the highly concentrated halo core diminishes 
until it disappears completely when xs ~ Xs.crit ~ 490xs,o 
(for T4 = 2). 



2.4 Evaporation time 

A rough, order-of-magnitude guide to the time-scale t^v for 
minihalo photoevaporation is the sound-crossing time for the 
characteristic size of the minihalo at the sound speed of the 
ionized gas, tsc = 2r-t/cs(10''K). The actual I-front propaga- 
tion and evaporative wind are quite complex, as we shall see, 
and depend upon additional properties like the halo mass 
and density, flux level and spectrum, and the time it takes 
the halo and background universe to evolve. However, this 
simple estimate will be a convenient standard of reference 
with which to compare the actual simulation results. 

For a minihalo of mass M at high redshift, we obtain 

t.c = 98Myr(M0^/^(^^) '^\l + z);,\ (24) 

using the adiabatic sound speed Cs(10'*K) = 

11.7(T4/m)^''^ = 15.2Ty\ms"\ where fi = 0.59 is 
the mean molecular weight for fully ionized gas of H and He 
if the abundance of He is v4(He) = 0.08 by number relative 
to H, and we have assumed T4 = 1. Of course, if tcv ^ tn, 
the Hubble time at the epoch of photoevaporation, then 
the underlying properties of the halo, infall, ionizing source, 
and background universe will evolve significantly during the 
process. The Hubble time at high redshift in the ACDM 
universe is 

^--3ii)-^^^^^^(w)"'^^ + ^)-"- 
Thus, the ratio of the two timescales is 

fl = 0.18HMrr^(^^^y\l + z)li\ (26) 

We see that this ratio is only weakly dependent on the back- 
ground cosmology, and suggests that the photoevaporation 
of a minihalo over the range of mass and epochs of interest 
will be completed in less than a Hubble time, if this estimate 
of tcv is correct. 
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Figure 4. Initial conditions for the minihalo photoevaporation 
simulations (TIS halo + self-similar cosmological infall): atomic 
number density (upper panel), and radial velocity (lower panel) 
profiles vs. radius for minihalo of mass My = 1 at (1 + 2)10 = 1. 



2.5 Ionizing photon consumption 

An important quantity characterizing the effect of minihalos 
on the process of reionization is the consumption of ioniz- 
ing photons which results from photoabsorption by minihalo 
atoms during their photoevaporation. We can parametrize 
this consumption by the photon-to-atom ratio, ^ = N-y/Na, 
i.e. how many photons per minihalo atom are required to 
photoevaporate the halo. Collapsed minihalos are signifi- 
cantly denser than the IGM, thus the recombination rate 
inside is significantly higher, leading to a much higher pho- 
ton consumption rate per unit time per atom, which could in 
turn have a substantial eff'ect on the progress and duration 
of reionization. 

We define the effective absorption cross-section at fre- 
quency of a minihalo at a given time t as 



o"cff,i/(i) = 



27rrfl — e 



(27) 



where = t^.h + t^,hoI + t^.HoII is the total bound- free 
opacity at impact parameter r and time t, including the 
minihalo gas which is leaving the halo in a wind, where 



(Tht,i{y)nidx. 



(28) 



where (Jbf,i and t-th.i are the bound-free cross-section and 
ionization threshold frequency and Ui is the number den- 
sity for species i. The number of photons absorbed by the 
minihalo gas per minihalo atom over time tev is then given 
by 



1 



(29) 
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Figure 5. Initial time-slice at t = Myrs (2 = 9 in ACDM 
universe), from top to bottom: (1) isocontours of atomic density, 
logarithmically spaced, in (r, x)— plane of cylindrical coordinates; 
(2) flow velocities; arrows are plotted with length proportional 
to gas velocity. An arrow of length equal to the spacing between 
arrows has velocity 5kms~^; the minimum velocities plotted are 
Ikms"^. Solid line indicates the boundary of the minihalo; cuts 
along the r = axis of: (3) gas number density n, (4) velocity , 
and (5) Isothermal Mach number Mj. Source is located outside 
the box, far to the left along the x-axis. 



where Fv(t) is the time dependent ionizing photon fiux per 
unit frequency. 

Alternatively, ^ can be calculated by the direct count 
of the number of recombinations experienced by each atom 
which was initially in the halo: 



1 + 



1 



dt I d\/(a^^nenHii + a'i^gnenHcii), (30) 



(2), 



where rie is the number density of electrons, nnii and riHcii 
are the number densities of H II and He II, and a\j and Q-jj], 
are the Case B recombination coefficients for H II and He II, 
respectively, and the volume integral is over all Lagrangian 
fluid elements initially inside the minihalo when the I-front 
flrst encountered it. We have neglected the recombinations 
of He HI to He II because these generally contribute diffuse 
flux which is absorbed on-the-spot by H and He. 

A naive estimate of ^ is obtained if we assume that the 
minihalo is optically-thin and instantaneously ionized but 
remains static at its initial density for a time t^c- Then, 
ignoring the contribution of He to recombinations, equa- 
tion 1I3OI I yields 
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Cint{nH)a„ 



(2) 



1 + 5tis 

jHaiman. Abel fc MadaiJl200ll) . where {uh) 
atom number density inside a halo, Cint = 
444^^ is the effective clumping factor for 
1 + <5tis = 130.6 is the average overdensity 
with respect to the cosmic mean background 
in reality, such a model is oversimplified, we 
an "efficiency" factor / which we will derive 
lations below. Thus we obtain 



(31) 

is the mean H 
{nli)/{n„}^ = 
the TIS, and 
of a TIS halo 
density. Since, 
also introduce 
from our simu- 



^ = 206fT--"*M^ 



0.15 



(32) 



According to this estimate, if f fx 1 as claimed by 
iHaiman. Abel fc Madaul ll200j) . based on their simulation 
of a uniformly-ionized minihalo with zero optical depth, the 
photoevaporation of minihalos can be an enormous sink of 
ionizing photons during the reionization epoch. By contrast, 
atoms in the IGM at closer to the mean baryon density 
have otherwise been estimated previously to consume only 
^ ~ 1 photons per atom during reioniz ation llGnedin|l2000t 
iMiralda-Escude. Haehnelt fc R,eesl200flh . It is crucial, there- 
fore, for us to perform the detailed simulations reported here 
in order to derive the efficiency factor / in equation 13211 
properly. 



3 THE CALCULATION 
3.1 Basic equations 

3.1.1 Gas dynamical conservation equations 

We solve the Eulerian conservation equations of fluid dy- 
namics for the vector U of the densities of conserved quan- 
tities in 2-D, cylindrical symmetry. 



(33) 



mass density p, momentum densities pv^ in the a;-direction 
and pVr in the radial direction, energy density E — p{v^ -\- 
Vr)/2 + e and number density of chemical species i in ion- 
ization stage z, ni,z. Here Vx and v,. are the velocity com- 
ponents, e = cvT is the specific internal energy, T is the 
temperature and cv is the isochoric specific heat. The Eu- 
lerian conservation equations can then be written as 



dt dx dr ' 



(34) 



where F(U) and G(U) are the fluxes of U in the x and r 
directions, respectively, given by: 

F(U) = {pvx,pv1 + p, pVxVr,Vx{E + p),ni,zVx), (35) 
and 

G(U) = {pVr,pVxVr,pvl + p, Vr {E + p) , Ui^^Vr) , (36) 

and S are the source terms; 



pVr pVxVr 



(V0). 



(V0). 




MM 



2 4 6 

t (Myr) 

Figure 6. R-type I-front speed. Position xj (in units of the box 
length Xbox = 5.64 Mpc; upper panel) and velocity V[ (in units 
of the analytical prediction as described in the text; lower panel) 
of an R-type ionization front propagating through the mean IGM 
at redshift z = 9. 



VriE+p) 



V- V<?!>-A + r, 



-Sr. 



(37) 



Here p is the pressure, (/> is the gravitational potential, A is 
the cooling rate, F is the heating rate, and 5*^,2 is the source 
function for the chemical species i in ionization stage z, as 
described in § ISTO We use the ideal gas equation of state 
for monatomic gas (7 = 5/3) 

p = -^ep. (38) 



3.1.2 Nonequilibrium chemical reaction network 

We take account of the nonequilibrium ionization balance of 
H, He, and a possible ad mixture of heavy el e ments C, N, O, 
Ne, an d S, as described in lRaga et alJ lITQQTli . lMellema et alJ 
il998r . and references therein. [For a description of the de- 
tails and checks of our microphysics, the reader is referred 
to Mellcma (1993^]. 

The source function for the nonequilibrium chemical re- 
action network is given by 

— UeUi^zCi^z + ni,z-10i,z-l — ni.zfjii.z, (39) 

where, as above, i labels the chemical elements, z labels the 
ionization state, Qi,z+i(T) is the recombination rate from 
stage 2 -I- 1 to 2, Ci,z_i is the coUisional ionization rate 
from 2 — 1 to 2, — {F^cTv.i.z) / {hv)dv is the pho- 

toionization rate from 2 to 2 -f 1, and fth.i is the photoion- 
ization threshold for species i. Charge exchange reactions 
are included, although we have not displayed them in equa- 
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Figure 7. Detailed structure of R-type ionization front in the mean IGM at 2 = 9 at time 0.5 Myr after the ionizing fiux turns on at 
the left-side boundary of the box. (a) (left) BB 5e4 case and (b) (right) QSO case. From top to bottom; (1) isocontours of pressure, 
logarithmically spaced, in (r, x)— plane of cylindrical coordinates; cuts along the x-axis (the I-front propagation direction) of (2) pressure, 
(3) temperature T, (4) H I (solid) and H II (dotted) fractions; (5) He I (solid). He II (dashed) and Ho III (dotted) fractions; (6) bound- free 
optical depths measured from a; = along the x-axis, for H I (solid), He I (dashed), and He II (dotted) at their respective ionization 
thresholds. Source is located 0.5 Mpc, to the left of the box. 



t=0.5 Myr, BB le5,F„= 1 ,(0.3,0.7,0.7),w/I-front,lGM,z = 9 
P 



3x1022 




Figure 8. Same as in Figure 171 but for BB le5 case. 



tion 1391 1. Recombination rates include radiative and dielec- 
tronic rates. 

The electron number density rie is given by rie = ?ih 11 + 
JT-Hc II + 2nHe III + nc ■ The contribution of one electron per 
carbon atom is added to ensure that the electron number 
density is always positive. Since the ionization threshold of 
carbon is below 13.6 eV and the IGM is transparent to these 
photons even before reionization, C is always at least singly 
ionized. The contribution of the other metals to the electron 
density is assumed negligible. 



3.1.3 Radiative cooling and photo-heating 

Our nonequilibrium cooling function includes collisional line 
excitation and ionization, recombination and free-free cool- 
ing due to H, C, N, O, Ne, a nd S, as described and tested 
in detail in lRaga et alJ I^Q^), to which we have added He 
cooling and the Compton cooling which results from inverse 
Compton scattering of the CMB photons off free electrons. 
The Compton cooling rate is given by 



Ao 



5.65 X 10 



■'^ne(l + z)*(r- 



''(40) 



TcMB)ergcm "s 

JShapiro fc Kandll987D . where Tcmb = 2.73(1 + z) is the 
CMB temperature at redshift z, Ue is the local density of 
free electrons, and T is the local gas temperature. 
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The photo-heating rate T (in units of ergcm~'^s~^) is 
the sum of the ionization heating terms for each of the 
species H I, He I, and He II: 



hv 



-0'hi,v,\ 



(41) 



is chosen so that the turn-around radius of the cosmologi- 
cal infall starts well inside the box. This ensures that the 
gas next to the simulation boundaries (excluding the axis of 
symmetry) is outflowing throughout the simulation, consis- 
tent with our transmissive numerical boundary conditions 
(see §EJ- 



3.1.^ Radiative transfer and the ionizing flux 

The energy flux of ionizing photons at a point with coordi- 
nates (r, x) for a source on the s-axis at position xo is given 
by 



F,{r,x) 



(42) 



Here Lv is the source luminosity at frequency and r^(r, a;) 
is the bound-free optical depth from the source to the point 
{r,x), given by 



Tv{r,x) = dvMlNm. + 0-i,,HclAfHcI + O-^.HoIiA'hcII, 



(43) 



where A^^.^ is the column density of species i at ionization 
stage 2 along the line from the source to the point (r,x). 



3.1.5 Gravity force 

We have modified the code to include the gravity field due 
to the halo of dark matter and gas, and the infalling matter, 
as discussed in detail in the Appendix. We neglected the 
small change in this gravity force over time which results 
from the gas dynamical evolution which moves the gas rela- 
tive to the dark matter. Thus the density field used for the 
gravity calculation is spherically-symmetric at all times and 
the gravitational acceleration is given by 



^2 



(44) 



where M(^ R) is the time- varying mass within the spher- 
ical radius R = {r^ + x^)^^^. This force is included in the 
Euler equations 13411 through the source term S as shown in 
equation 1371 1. 



3.2 Initial conditions 

3.2.1 Minihalo and cosmological infall 

As discussed above, the halo in our illustrative simula- 
tions has radius rt — 0.76 kpc, and virial temperature 
Tvir = 4000 K, corresponding to dark-matter velocity dis- 
persion av ~ 5.2kms~^. The initial number density and 
velocity profiles of the halo and the infall are shown in Fig- 
ure 3] For comparison with our simulation results at later 
times, we show the initial conditions on the computational 
grid in Figure]^ (top to bottom): (1) the density contours 
(logarithmically-spaced), (2) the initial velocity field, and 
cuts along the axis of: (3) gas number density n, (4) veloc- 
ity Vx, and (5) Isothermal Mach number Mi = v/csj, where 
Csj is the isothermal sound speed. The computational box 



3.2.2 Source spectra and evolution 

The source is on the axis of symmetry, far away (thus the 
rays are close to parallel, although the code actually takes 
the true angles properly into account) to the left of the box. 
We consider three possible spectra as described in §0 two 
stellar cases, black-body spectra with effective temperatures 
T = 50,000 K ("BB 5e4") for Pop. II stars and T = 10^ 
K ("BB le5") for Pop. HI stars, and a power-law QSO- 
like energy spectrum with index ct = —1.8. As the universe 
expands, the proper distance between the source and the 
minihalo grows, as xoit) oc a{t), where a{t) is the cosmic 
scale factor, so the incident flux level decreases oc Xq^ oc 
a{t)-\ 



3.2.3 R-type I-front encounters the simulation volume at 
t = 

The I-front has a finite width of ~ (10 — 20) mean free paths. 
For the intergalactic weak, R-type I-front which encounters 
the minihalo, this mean free path in the IGM is a few kpc, 
which is similar to our simulation box size. Therefore, it 
would be inaccurate to start the simulations by assuming 
a sharp I-front at the left box boundary. Instead, we start 
our simulations shortly before the I-front arrival at the left 
boundary, as follows. We assume the the gas outside the box 
is uniform, at the mean IGM density, for which case there 
is an exact solution for the weak, R-type I-front propaga- 
tion. We then use the frequency-dependent optical depths 
given by the current position of the I-front obtained from 
the analytical solution to attenuate the spectrum of the pho- 
toionizirig source which enters the computational box. While 
not exact (e.g. due to infall which slightly increases the lo- 
cal density around the halo) this approach closely imitates 
the gradual rise in the photoionizing radiation fiux arriving 
from the source as the I-front approaches and the hardening 
of the spectrum due to the deeper penetration of higher- 
energy photons. As an alternative, it is also possible to take 
a much larger simulation box, which would allow the I-front 
to relax to its correct structure before arriving at the halo. 
However, such an approach would either degrade the res- 
olution unacceptably, or else the simulation would become 
prohibitively (and needlessly) more expensive. 



4 NUMERICAL METHOD AND TESTS 
4.1 The method 

We use a 2-D axisymmetric Adaptive Mesh Refinement 
(AMR) code (called CORAL) described in detafl in 
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Figure 9. Ionization structure of metals in R-type I-front propagating in the mean IGM at z = 9, 1.3 Myr after front enters the box on 
the left-hand-side. C, N, and O ionic fractions along symmetry axis: (a) (left) BB 5e4 case; (b) (right) QSO case. 
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Figure 10. Same as in Figure|5] but for BB le5 case. 



iRaga et alJ (Il995l) and iMellema et al.1 lll998h . The code has 
been extensively tested and applied to a range of problems 
including simulations of Herbig-Haro jets and the photoe- 
vaporation of uniform density interstellar clouds near plan- 



etary nebulae. Here we will discuss the numerical details of 
the code only briefly, concentrating on the new elements in 
the code introduced by us, and we refer the reader to these 
earlier papers and references below for a more in-depth dis- 
cussion. 

The numerical scheme for solving the Euler equations is 
an adaptive grid i mplementation of the van Leer Flux Vector 
Splitting method (EnLee3[l983), improved to second order 
accuracy by use of linear gradients within cells as described 
in Arthur ( 1991). The refinement and de-refinement criteria 
are based on the gradients of all code variables. When the 
gradient of any variable is larger than a pre-defined value 
the cell is refined, while when the criterion for refinement 
is not met the cell is de-refined. We have modified the code 
to include the gravitational acceleration in equation 1441 , as 
discussed in detail in the Appendix. We impose refiective 
boundary conditions on the axis and transmissive boundary 
conditions on the other three boundaries of the box. 

The code timestep is set in general by the minimum 
of the timesteps determined by the Courant condition, the 
ionization time scale and gravity. 

At = min(eCourAfCour, EionAti 

Oil) tgrav Afgrav). (45) 

The Courant time in equation 14511 is 

AfCour = Ax/cs , (46) 

where Aa; is the cell width. The ionization time scale here 
is that for ionizing hydrogen. 



Ation — ^io 



' (dnm/dt) 



nnifpHi — nt,,,a 



(2) I 



(47) 
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where 4>hi is the photoionization rate for hydrogen, tihi and 
riHii are the number densities of neutral and ionized hydro- 

(2) 

gen, ajj is the Case B recombination coefficient for hydro- 
gen, and eion is a constant to be determined. The gravity 
time is just 

xl/2 



At„ 



(Ax/ag 



(48) 



The constants ei are small numbers whose values are chosen 
to ensure accuracy and stability while minimizing the total 
number of timesteps per simulation. In all cases considered 
here, the timestep is set by the Courant condition, except 
when the fast, R-type I-front propagates through the com- 
putational box. In the latter case, the ionization time scale is 
smaller than the Courant time, and the value of eion is chosen 
by experimentation to optimize the accuracy and efficiency 
of the test problem described in S 14.31 For the Courant con- 
dition, the van Leer Flux Vector Splitting scheme is stable 
when ecour ^ [2y + M(3 - 7)]/(7 -|- 3), where 7 = 5/3 is 
the adiabatic index and M is the Mach number. For M = 
we obtain the condition ecour ^5/7 for the scheme to be 
stable for any Mach number M. For the current simulations 
we utilized the conservative value ecour = 0.4. 

The microphysical processes - chemical reactions, radia- 
tive processes, transfer of radiation, heating and cooling - 
are implemented though the standard approach of operator- 
splitting (i.e. solved each time-step, side-by-side with the 
hydrodynamics and coupled to it through the energy equa- 
tion). The energy and chemical rate equations are solved 
semi-implicitly. We follow the nonequilibrium evolution of 
the ionic species of H, He, C IT VI, N TVI, O TVI, Ne 
TVI, and S II- VI. The species C I and S I are assumed 
fully ionized since their ionization thresholds are below 
the ionization threshold of hydrogen and the gas is largely 
optically-thin t o suc h low energy photons [for details see 
[MgUc ma et al] lll998h and references therein]. The method 
for solving the nonequilibr i um ch emistry equations is from 
ISchmidt-Voigt fc KoePD enl (Il987f) and is imp lemented as de- 
scribe d and tested in lRagaetalT il997ft and lMellema et all 
l|l998 h. We use Case B recombination coefficients and the 
corresponding cooling rates, as is appropriate for the range 
of densities and temperatures in these simulations. 

Radiative transfer of the ionizing photons is treated ex- 
phcitly by taking into account the bound-free opacity of H 
and He in the photoionization rates and heating rates as ex- 
plained in detail in Mell ema et al., (1_998). In order to speed 
up the radiative transfer calc ulation, the optical depths are 
approximated as described in lTenorio-Taele et alJ il983h . as 
follows: 



{r„,A for z^th.Hi ^ «^ ^ fth.Hci 
T^,B for Z^th.HcI ^ ^ i^th,HcII 
r^,c for Vth.HcU ^ 1^, 



where 
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Figure 11. Evolution of the I-front from R-type to D-type inside 
the minihalo for BB 5e4 case, (upper panel) Position xi (relative 
to minihalo center, in units of the minihalo radius rt at t = 0) and 
(middle panel) velocity vi of the I-front as it travels toward and 
across the minihalo. The positions of the halo boundary (long- 
dashed line) and centre (short-dashed line) are also indicated, 
(lower panel) Fraction of mass Mi^it, the mass which is initially 
inside the minihalo when the intergalactic I-front overtakes it, 
which remains neutral versus time t. The arrows in the top and 
bottom panels mark the analytically-calculated width of the ISL 
for zero impact parameter and the mass fraction inside the ISL, 
respectively, while the arrow on middle panel marks the moment 
when I-front becomes R-critical, defined by the condition tij = 
2Csj 2i the onset of transition from R-type to D-type, where / 
is isothermal sound speed of ionized, post-front gas. 



-h(1.81)''VHel] 



(51) 



and 



+4 r^.Hoii 



(52) 



Tv,B = [(0.63)^''^r„,Hi -I- 



Here Ti^^i^z ~ cr,y,th,i,zA'i are the optical depths at the respec- 
tive ionization thresholds of H I, He I and He II. This allows 
us to precompute the integrals over frequency involved in 
calculating photoionization and photoheating rates for each 
assumed source spectrum to make an extensive look-up table 
for these rates as functions of the threshold optical depths, 
for use in the rapid computation of these rates during each 
run by interpolation between entries in the look-up table. In 
order to speed up the calculations further, the optical depths 
are normally re-calculated once every 5 time steps. 
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Figure 12. Weak, R-type phase of I-front in IGM just about to overtake a 10^ Mq minihalo at z = 9, t = 0.2 Myr after the I-front 
entered the box at left-hand side (ionizing source is located far to the left of computational box along the a;-a3cis): (a) (left) BB 5e4 case 
and (b) (right) QSO case. Panels show the same quantities as in Figure |7| The current position of the I-front (50 % ionization level of 
hydrogen) is indicated on the uppermost panel with the dashed line. 



4.2 Simulation parameters 

Our box size in {r,x) dimensions is 3.45 kpcx6.9 kpc, except 
for the uniform IGM simulations in § 14.31 and § 15.11 where 
we adopted larger boxes, as described there. The resolution 
at the finest grid level (i.e. fully-refined) is 1024x2048 cells, 
again with the exception of the simulations in ii l4.3l and i) l5.1l 
where such high resolution was not required, but a longer 
propagation time was necessary, so we used resolutions of 
64 X 8192 and 512x1024, respectively. 

The gas consists of hydrogen and helium in primordial 
abundance (24.2 % He by mass), with a small trace of metals 
(C, N, O, Ne, and S) at 10"^ Zq. The initial ionizing fiux 
we adopt is = 1, except in 14.31 where we use Fo = A. 
Each photoevaporation simulation took from several days to 
1 week on a 2 GHz Athlon processor. 

4.3 Test problem: R-type I-front in uniform IGM 

In order to track properly the supersonic, weak, R-type I- 
fronts which sweep into and around our minihalos, we de- 
fined an "ionization time-step" related to the time-scale for 
ionizing hydrogen, as described above in ^ 14. II We tested the 
scheme and adjusted the value of eion by performing a simu- 
lation of an ionization front propagating through a uniform 
density IGM in a long and thin (44 kpc x 5.64 Mpc), quasi- 
ID simulation box at resolution 64 x 8192. In this case, the 
gas is pure hydrogen, so there is a well-k nown exact analyt- 
ical solution for the I-front propagation iShapiro fc Girou^ 



The density of the gas is equal to the mean cosmic 
baryon density at z = 9. The source is located on the sym- 
metry axis, 0.5 Mpc to the left of the computational box, 
with a flux as measured at the left boundary of the box equal 
to Fo = 4. At time t = Q, the box is neutral. We identify the 
location of the I-front as the position at which the ionized 
fraction is x = 0.5. The result is plotted in Figure 6. By 
experiment, we found that eion = 0.05 is sufficient to obtain 
the correct velocity for the R-type I-front to better than 5 
%, as shown in Figure |S] Using larger eion can lead to an 
underestimate of the velocity by a significant amount. 



5 RESULTS 

5.1 Structure of global I-front in the IGM during 
reionization: the weak, R-type phase 

To study the detailed temperature and ionization structure 
of the global I-front during reionization in its weak, R-type 
phase and the dependence of these on the source spectrum, 
we performed a set of three simulations of an I-front in the 
uniform IGM at z = 9, by "zooming in" with a smaller box 
(34.5 kpc X 69 kpc) and better length resolution (512x1024, 
fully-refined) than were used in the test problem in the pre- 
vious section. The source in each case produced a flux Fo = 1 
at the left boundary of the simulation box, with a different 
spectrum for each simulation. The results are shown in Fig- 
ures ITI llOl As expected, for the cases with harder spectra 



© 0000 RAS, MNRAS 000, 000-000 



Photoevaporation of Cosmological Minihalos during Reionization 17 



t = 0,gO Myr, DP le5 ( 1 O^M^.O. l),(0.3.0.7,0.7).w/I-front 
P 






a: 




Vv^-^r^-,-, 



: logio(T„) ' 









X (cm) 



2x1022 



Figure 13. Same as Figure [T^ but for BB le5 case. 



(i.e. QSO, BB le5), the I-front is broader than in the case 
BB 5e4, due to the more deeply penetrating hard photons 
of the former cases. There is significant pre-heating and pre- 
ionization in the former cases, ahead of the point inside the 
I-front at which the H neutral fraction is 0.5. In all cases, the 
ionization of He is predominantly from He I on the neutral 
side to He HI on the ionized side, with only a small fraction 
of He II. The position of the hydrogen I-front (sh = 0.5) is 
nearly coincident with the corresponding point for helium 
(a^Hc = 0.5) for the softer BB 5e4 spectrum, while for the 
harder QSO and BB le5 spectra the He I-front is more ad- 
vanced by ~ 7 kpc. Another important difference to note is 
that the post-front temperature is ~ 2 x 10^ A" in the BB 
5e4 case, but ~ 3 x 10* A" in the QSO case, and even higher, 
~ 4 X 10* if in the BB le5 case. The ionization states of the 
trace of heavier elements also differ significantly for differ- 
ent source spectra, as shown in Figures 1^ and [THl The BB 
5e4 spectrum ionizes C, N, and O mostly to C HI, N HI, 
and O HI, respectively, with only small fractions, of order 
10 %, of C IV, N IV, and O IV, while the BB le5 spectrum 
ionizes these species mostly up to ionization stage IV, with 
a notable fraction (> 10%) of ionization stage V, and tiny 
fractions of O VI and N VI. Finally, the hard QSO spectrum 
is able also to produce fractions as high as a few per cent 
for O VI and N VI, while the bulk of the metals reside in 
the ionization stages IV and V. 



5.2 I-front encounters the minihalo: weak, R-type 
phase inside minihalo 

As the weak, R-type I-front propagates towards the mini- 
halo, it encounters a steeply rising density, first the infall 
region and then the virialized region of the minihalo, itself, 
and its velocity drops precipitously. Figure 1111 shows the 
evolution of the I-front position xi (upper panel), its veloc- 
ity vi (middle panel) and the mass fraction M/Mi-nn of the 
gas initially within the original hydrostatic sphere (bottom 
panel) which remains neutral over time in the BB 5e4 case. 
The initial speed of the I-front as it enters the box is about 
10,000 kms~^, dropping to ~ 1000 kms~^ by the time the 
front crosses the virial radius of the halo. As a result, the 
I-front is initially a weak, R-type front even after it enters 
the minihalo. It takes about 5 Myr for the front to slow to 
the R-critical front speed vr = 2cs,2 ~ 20kms~^ (where 
the isothermal, post-front sound speed Cs,2 — llkms"^) 
(marked on the plot by an arrow) after which a compressive 
shock must form to lead the I-front and further decelerate it, 
so as to transform it to a D-type front. On the top panel, we 
have indicated (by another arrow) the position of the inverse 
Stromgren surface at zero impact parameter, as calculated in 
§12131 The two arrows almost coincide, confirming our expec- 
tations of approximate correspondence between the ISL (in 
the static analytical calculation) and the surface on which 
the I-front is R-critical, close to the moment of conversion 
of the I-front to D-type (also indicated in Fig.|5J. This cor- 
respondence is not perfect, however, due to the significant 
approximations made in the ISL calculation. For example, 
due to the centrally-concentrated halo density profile the I- 
front transition to D-type occurs in fact at different times for 
each impact parameter, leading to a continuously evolving 
I-front shape different from the shapes shown in Figure |21 

Figures 1121 and 1131 show the time-slice at 0.2 Myr, 
shortly before the weak, R-type I-front enters the minihalo. 
The instantaneous position of the I-front is indicated on 
the top panel which shows the contours of pressure. The 
I-front moves faster further away from the halo due to the 
lower density in the infall profile and has already signifi- 
cantly slowed down closer to the axis. The harder photons 
clearly penetrate significantly deeper in the QSO and BB 
le5 cases and the halo shadow is already clearly outlined at 
that time by the pressure contours, even though the front 
itself has not passed the halo yet. These same hard photons 
also start to heat the halo on the source side, while deeper 
into the halo and in the shadow region behind it virtually 
no photons can penetrate due to the extremely high optical 
depth, r ~ 10* (bottom panels). While the shadow is mostly 
neutral in all cases, the QSO case produces a 10^* fraction 
of He II and He HI there. 

In Figures I14I17I we show the structure of the flow at 
t = 2.5 Myr, when the I-front is inside the minihalo but still 
a weak, R-type front, before much hydrodynamical back- 
reaction has begun. Figure 1141 and 1151 show that the gas 
which was initially infalling on the source side next to the 
halo is already reversing its velocity, however, and starting 
to form a shock which will sweep the IGM outward. On the 
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Figure 14. Weak, R-type phase of I-front inside the minihalo. Same as Figure 5, but for later time-sUce, t = 2.5 Myr after I-front enters 
the box: (a) (left) BB 5e4 case and (b) (right) QSO case. A velocity arrow of length equal to the spacing between arrows has velocity 
lOkms"^; minimum velocities plotted are 3kms~^. Solid line shows current extent of gas initially inside minihalo at z = 9. Dashed line 
indicates the current position of the I-front (50% H-ionization contour). 
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Figure 15. Same as in Figure [T^ but for BB le5 case. 



other hand, on the shadow side of the halo, the gas infall 
continues uninterrupted. Also evident is some squeezing of 
the gas in the shadow due to the much higher pressure of 



the ionized gas outside it after the passage of the fast R-type 
I-front there. 



5.3 I-front trapping: R-critical phase and the 
transition from R-type to D-type 

The I-front finally slows to become an R-critical front by 
5 Myr. We show the results for this time-slice for all three 
source spectra in Figures I18H22I To illustrate clearly that 
this is the R-critical phase, we focus on the Pop II stel- 
lar case BB5e4 in Figure 1181 and "zoom in" on the region 
around the I-front along the a;-axis in Figure 1191 All of 
the important properties of an R-critical I-front are evi- 
dent in Figure 1191 The I-front velocity at this epoch is 
?;/ « 16 km s~^ which is not far from the estimated R-critical 
value of vr ~ 20 kms~^, calculated based upon the immedi- 
ate post-front temperature of T2 ~ 10, 000 K, the pre-front 
Ti = 4,000K (i.e. the neutral side is undisturbed mini- 
halo gas at this virial T), ^1 = 1.11 (i.e. 90% neutral), and 
^2 = 0.63 (i.e. 90% ionized H), which yield the post-front 
and pre-front isothermal sound speeds Cs,/,2 = llkms"^ 
and Cs,/,i = 5.3kms~^, respectively. The numerical results 
show the predicted jump in density by a factor of close 
to two across the R-critical I-front, required by the I-front 
jump conditions. In addition, in the lab frame in which the 
neutral minihalo gas on the pre-front side is at rest (i.e. 
the frame in which we have plotted velocity and isothermal 
Mach number in Figures il8l and [T^ . the ionized gas just 
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Figure 16. Weak, R-type phase of I-front inside the minihalo. Same as 

Figure El but for t = 2.5 Myr. 
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Figure 17. Same as in Figure [T?)l but for BB le5 case. 



behind the front should be moving toward the neutral side 
at the isothermal sound speed of the ionized post-front gas, 
or about llkms"^, which our plots confirm (i.e. note that 
the isothermal Mach number Mj = 1 at that point in lower 
panel of Fig. 1191 . It is interesting to note that, while the gas 
on the ionized side of the I-front is generally optically thin 



during most of the evolution, the optical depth is not neg- 
ligible a,t t — 5 Myr. According to Figure \W\ the H optical 
depth across the ionized, post-front layer is of order a few, 
so the speed of the I-front is lower than would be calculated 
using the unattenuated flux in the I-front continuity jump 
condition. 

This R-critical phase is the beginning of the transition 
to D-type, which should occur shortly, thereafter, when a 
shock wave forms to compress the gas ahead of the front, to 
enable the front velocity to drop to the D-crit value. In this 
case, using the sound speeds at 5 Myr reported above, the D- 
critical velocity is vd ~ c^j i/(2csj,2) — 1.35 kms^^ ^ vr. 
As we can see from the plot of vi in Figure ITTI however, the 
transition from R-type to D-type is quite extended in time, 
lasting lO's of Myr. During that same time, the I-front ad- 
vances into denser and denser neutral minihalo gas, which 
also slows the front quite apart from any dynamical com- 
pression which might lead the front. In the meantime, from 
the R-critical phase onward, the hydrodynamical response 
to the I-front becomes more and more dramatic. 



5.4 The structure of the photoevaporative flow 
during the D-type I-front phase 

When the I-front slows to become a D-type front, the hydro- 
dynamical response of the gas catches up with it and leads 
to the complete photoevaporation of the minihalo. This pho- 
toevaporative flow exhibits many generic features which are 
independent of the spectrum of the ionizing source. The side 
facing the source expels a supersonic wind backwards to- 
wards the source. There is a wind shock which thermalizes 
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Figure 18. R-critical phase of I-front just before transition from R-type to D-type. (a) (left) Panels show same quantities as in Fig. Ill 
but for t = 5 Myr and BB5e4 case; (b) (right) Same quantities as in Fig. 1121 but for the t = 5 Myr and BB5e4 case. 
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Figure 19. Same as Fig. 1181 but "zoom-in" to enlarge the view of R-critical I-front along the x-axis. 



this wind outflow and separates the unshocked supersonic 
wind close to the ionized side of the I-front from the shell of 
shocked wind further away from the I-front on this side. The 
shocked wind gas acts as a piston which sweeps up the pho- 
toionized IGM outside the minihalo and drives a shock into 



the IGM, reversing its infall velocity. The wind grows more 
isotropic with time as the remaining neutral halo material 
is photoevaporated. Since this gas was initially bound to a 
dark halo with av < lOkms"^, photoevaporation proceeds 
unimpeded by gravity. In Figures 1^^261 we show the struc- 



© 0000 RAS, MNRAS 000, 000-000 



Photoevaporation of Cosmological Minihalos during Reionization 21 



1=5 Myr. QSO. (10'Mg,9.1), (0.3.0.7,0.7). (10. 204Hx1024 t=5 Myr. BBleS. (10'Mg,9.1), (0.3.0.7,0.7). (10.3). 204Hx1024 




X (cm) X (cm) 



Figure 20. R-critical phase of I-front just before transition from R-type to D-type. Panels show same quantities as in Fig. 1141 but for 
5 Myr; (a) (left) QSO case; (b) BBleS case. 
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Figure 21. R-critical phase of I-front just before transition from R-type to D-type. Panels show same quantities as in Fig. 1121 except 
for t = 5Myr: (a) (left) QSO case; (b) BBle5 case. 
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Figure 22. D-type phase of I-front: Photoevaporative Flow 60 Myr after I-front overtakes the minihalo: (a) (left) BB 5e4 case and (b) 
(right) QSO case. Same quantities plotted as in Fig. 1141 except t = 60 Myr. 
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Figure 23. (a) (left) Same as Figure l22l but for BB le5 case, (b) (right) Same time-slice as Figure l22l but for optically-thin BB 5e4 
case. 
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Figure 24. D-type phase of I-front: Photoevaporative Flow 60 Myr after I-front overtakes the minihalo: (a) (left) BB 5e4 case and (b) 
(right) QSO case. Panels show the same quantities as in Figure [T^ Key features of the flow are indicated by the numbers which label 
them on the temperature plots: 1 = IGM shock; 2 = contact discontinuity which separates shocked halo wind (between 2 and 3) from 
swept-up IGM (between 1 and 2); 3 = wind shock; between 3 and 4 = supersonic wind; 4 = I-front; 5 = boundary of gas initially inside 
minihalo at 2 = 9; 6 = shock in shadow region caused by compression of shadow gas by shock-heated gas outside shadow. 




Figure 25. (a) (left) Same as Figure 0^ but for BB le5 case, (b) (right) Same as Figure 0^ but for optically-thin BB 5e4 case. 
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Figure 26. Snapshots of the temperature at times t = (top left), 0.2 (top right), 2.5 (middle left), 10 (middle right), 60 (bottom 
left), and 150 (bottom right) Myr in the BB5e4 case. Shaded isocontours are logarithmic in temperature. Colors indicate the values of 
logj^Q(T), as labelled on the color bar. 



ture of the flow at a time 60 Myr after the global I-front first 
overtakes the minihalo. In Figure 1^ and the left panels of 
Figure 123 we show the density contours of the gas, the cur- 
rent position of the ionization front, the velocity field, the 
current extent of the original halo material, and profiles of 
gas number density, velocity and Mach number along the x- 
axis. In all three cases, the typical speed of the photoionized 
outflow is supersonic and roughly proportional to the sound 



speed of the ionized gas. Hence, the outflow is fastest in the 
BB le5 case, due to the higher photoionization temperature 
reached in this case. 

For comparison, we have also performed a simulation 
of this problem in the optically-thin approximation (i.e. 
zero optical depth), to demonstrate the inadequacy of such 
an approximation. In Figure 1231 (right panels), we plot 
the results for this optically-thin photoevaporation simula- 
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Figure 27. Evolution of I-front position and velocity and of the of neutral gas content of photoevaporating minihalo: (a) (left) BB 5e4 
case and (b) (right) QSO case, (upper panels) Position xj (in units of the minihalo radius rt at t = 0) and (middle panels) velocity vj 
of an ionization front propagating toward and through the minihalo. The positions of the boundary (long-dashed line) and centre of the 
halo (short-dashed line) are also indicated, (lower panels) Fraction of mass Mi^it, the mass which is initially inside the minihalo when 
the intergalactic I-front overtakes it, which remains neutral versus time t. 



tion for exactly the same halo and source parameters as 
the BB 5e4 spectrum case shown in Figure 1221 (left pan- 
els), for the same time-slice. In the optically-thin case, the 
halo is suddenly and uniformly ionized and photoheated to 
T ~ 20, 000 — 35, 000 K, well above its virial temperature. 
Thereafter, the large pressure gradient causes the gas to ex- 
pand isotropically in all directions unimpeded by gravity, 
producing a completely different flow structure. 

For the more realistic simulations which included ra- 
diative transfer. Figures 1241 and 1251 show the pressure con- 
tours, and the profiles along the symmetry axis for pressure, 
temperature, H I-II fractions. He I-III fractions, and opti- 
cal depths at the ionization thresholds of H I, He I, and He 
H, for the photoevaporative flow 60 Myrs after the global 
I-front first overtakes the minihalo. Key features of the fiow 
are indicated by the labels on the temperature panels. For 
comparison, we also show the very different results for these 
quantities in the optically-thin simulations (Figure 1^ right 
panels) . 

5.5 Evolution of the temperature structure 

The different stages of the photoevaporation process are il- 
lustrated in Figure l^ bv shaded isocontours of temperature 
for a sequence of time-slices of the BB 5e4 case. The pan- 
els show: (a) (time t — Myr), the initial condition of the 
isothermal, virialized halo surrounded by a cold, infalling, 
IGM; (b) (t = 0.2 Myr), the fast, R-type I-front sweeping 



through the computational box is just about to encounter 
the minihalo along the x-axis, and is starting to photoheat 
the halo on the source side; (c) {t = 2.5 Myr), when the 
weak, R-type I-front is inside the minihalo before it slows 
to R-critical and begins to transform to D-type, a distinct 
shadow is apparent behind the neutral portion of the halo, 
while the rest of the gas is heated to T > 10'' K, and the pho- 
toevaporation process begins; (d) {t = 10 Myr), the R-type 
to D-type transition phase, when the I-front slowly advances 
along the axis, the shadow region behind the shielded, neu- 
tral minihalo core is reduced, and the IGM shock and the 
shocked wind are clearly visible; (e) {t = 60 Myr), strong 
D-type phase, when the shadow is being ablated and the re- 
maining gas heated, by shocks reflected off the a;-axis behind 
the halo. The key features of the flow indicated in Figure 1^ 
are all clearly seen; (f) {t = tov ~ 150 Myr, the evaporation 
time), when the gas is almost completely evaporated from 
the halo and resides in a much larger, low-density region 
which is cooling adiabatically as it expands, leaving behind 
a dark-matter halo almost completely devoid of gas. Movies 
of these simulation results, including the time evolution of 
additional quantities like pressure, density, H I and He I 
fractions, and the results for the cases with different source 
spectra are available at http://galileo.as.utexas.edu. 
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Figure 28. Same as Figure l27l but for BB le5 case. 
5.6 I-front evolution 

The evolution of the position xj and velocity vj of the I- 
front, and of the remaining neutral mass fraction M/Minit 
are shown for all three source spectra cases in Figures l?7l and 
1281 for comparison. The weak, R-type I-front slows down to 
the R-critical velocity about 5 Myr after it reaches the mini- 
halo, after which a shock must form ahead of the I-front to 
compress the gas and slow it down to D-critical velocity or 
below, as discussed in detail in S 15.21 The further evolution 
is slower and more gradual. Eventually, the I-front speed 
drops below the D-critical front speed vd ~ 1km s~ . For 
the rest of the evaporation process, the velocity of the I-front 
is < lkms~^, i.e. sub-critical. At 60 Myr, for example, the 
I-front velocity has dropped to about vi ~ 0.3kms~^, while 
the D-critical velocity is still vi ~ Ikms"^. This behaviour 
is different from the usual approximation made in analyti- 
cal calculations of photoevaporation, in which the I-front is 
assumed to be D-critical. 

What kind of D-type front is it during the D-type 
phase? Weak, D-type fronts move subsonically with respect 
to both the neutral and ionized gas. Strong D-type fronts 
move subsonically into the neutral gas but supersonically 
with respect to the ionized gas. In the lab frame (i.e. rest 
frame of the neutral gas before ionization), as the D-type 
front advances into the neutral side (i.e. away from the 
source), the ionized gas behind it always moves toward the 
source (i.e. opposite to direction of I-front). This ionized 
gas motion can be either subsonic (for weak fronts), or al- 
most sonic (for D-critical fronts), but it is supersonic only 
for strong D-type fronts. If we apply this description to our 
60 Myr time-slice for the BB5e4 case and look both at the I- 



front velocity evolution plot in Figure ITfl and the cuts along 
the i-axis in Figure VFR we see that the I-front is subcritical 
at that point. From the Mach number plot at points just be- 
hind the I-front on the immediate post-front (ionized) side, 
the lab frame velocity is supersonic toward the source, which 
indicates that the front is a strong, D-type. This is consis- 
tent with the density cut along the i-axis which shows that 
density drops by a very large amount from the neutral to the 
ionized side, as it should for both weak and strong, D-type 
fronts, but the strong type has an even bigger drop in den- 
sity, as required to explain the numerical results. In short, 
the I-front at 60 Myr is a subcritical, strong, D-type. 

5.7 Evaporation times 

The neutral mass fraction of the original minihalo gas 
quickly drops to ~ 60 per cent during the initial slowing- 
down phase of the I-front evolution, and subsequently de- 
clines much more gradually as the minihalo photoevapo- 
rates. We define the photoevaporation time tov as the time 
when only 0.1 per cent of the mass remains neutral (i.e. 
when M/Minit = 10"^). We see from Figures 1771 and 1^ 
that, once the neutral mass fraction reaches ~ 1 per cent, 
it drops fairly precipitously, indicating that the value of tev 
is not sensitive to our particular choice of the final value 
for Af/Afinit used to define this evaporation time. We obtain 
tcv ~ 150, 125 and 100 Myrs for the BB 5e4, BB le5 and 
QSO cases, respectively. Compared to the BB 5e4 case, the 
evaporation time in the BB le5 case is shorter due to the 
higher outflow velocity in this case (see Fig. I23II . caused by 
the higher temperature in the photoionized region. In the 
QSO case, on the other hand, the outflow velocity is similar 
to the one in the BB 5e4 case, but the photoevaporation pro- 
cess in this case is accelerated by the significant pre-heating 
due to the hard photons present in this spectrum, which 
ultimately leads to an even shorter tcv. These values of fov 
are all within a factor of two of the sound crossing time for 
an ionized gas sphere of the same diameter as our minihalo 
reported in §2.4. However, we caution that this result is not 
a proof that tsc is a good estimator of t^v In fact, results for 
a wide range of halo and source parameters which we shall 
present in a companion paper demonstrate that tcv departs 
significantly from t^c, in general. 

5.8 Ionizing photon consumption 

5.8.1 Minihalo effective cross-section 

The gradual decay of the opaque cross section of the mini- 
halo as seen by the source [as defined in eq. ll?7H is illus- 
trated by Figure l^ deft panels) for photons at the ionization 
thresholds of H I, He I, and He II. For the BB5e4 case, the 
cross section at the H I threshold drops to 10% of its original 
value after O.itev = 60 Myr, while the cross section at the 
He I threshold drops somewhat faster. The cross section of 
the minihalo at the He II threshold, however, initially rises 
as the partially ionized gas expands, and only later, around 



© 0000 RAS, MNRAS 000, 000-000 



Photoevaporation of Cosmological Minihalos during Reionization 27 



All spectra, ( ICM^.g, 1), (0,3,0.7,0.7) 




1 = 5,25.50,100 & 150 Myr, BB 5e4, (10'Mj,,9.1), (0.3.0.7,0.7) 

— ^ < \ > > \ r 




3x10^1 6x1021 
r (cm) 



9x1021 



Figure 29. (a) (left) Fraction of minihalo initial geometric cross section nr^ (rt = 0.76 kpc) which is opaque to source photons that can 
ionize H I (top panel), He I (middle panel), or He II (bottom panel) versus time (in units of tev), for the BB 5e4 (filled squares), QSO 
(filled triangles) and BB le5 (open triangles) cases, (t^v = 150,100, and 125 Myr, respectively), (b) (right) Transmitted flux Ii^/I^fi at 
the ionization thresholds of H (top), He I (middle), and He II (bottom) vs. impact parameter (or radius in axisymmetry) r at times t = 5 
(solid), 25 (dotted), 50 (short-dash), 100 (long-dash), and 150 (dot-dash) Myr for the BB 5e4 case. 



OAtev = 60 Myr, does it start to decline slowly, as He II 
starts to become ionized to He HI. The evolution of the H 
I effective cross-section is similar for the two stellar spec- 
tra, but decreases somewhat more slowly in the QSO case, 
while the evolution of the He H cross-section is similar for all 
three spectra. Finally, the evolution of the He I cross section 
is markedly different in all three cases, due to the different 
photon energy distribution of the three spectra. 

In Figure 1291 we show the transmitted fluxes at the 
ionization thresholds of H I, He I, and He II vs. impact pa- 
rameter r at times t = 5, 25, 50, 100, and 150 Myr for the 
BB 5e4 case. We see that the simulation box is optically- 
thin to the photons with energies below the He II ionization 
threshold along most lines of sight along the x-direction out- 
side the minihalo; only the minihalo is significantly opaque, 
once the global I-front has swept past the minihalo. The 
minihalo's cross-section declines quickly with time once the 
I-front enters it and is somewhat smaller for He I ionizing 
photons than for H I ionizing photons. The simulation box 
is never completely optically-thin to He II ionizing photons 
and has rather complex time behaviour as the wind expands 
and more He II is photoionized to He HI. 

5.8.2 Number of ionizing photons per atom consumed 
during photoevaporation 

Based on our simulation data, we have calculated the num- 
ber of ionizing photons per atom, ^, required to evapo- 



rate a minihalo both by using the effective cross-section 
method of equation 1291 1 and by counting the total number 
of recombinations [as in eq. (1301 1. Using the first method, 
we obtain ^ = (4.51,5.41,3.43) in the BB 5e4, QSO, and 
BB le5 cases, respectively, while the second method yields 
^ — (5.1, 5.0, 3.3). Since the two methods give similar results, 
we can, henceforth, rely upon either one. For the "efliciency" 
factor / in equation 1321 (assuming T4 = 1), these results 
imply / = (0.029, 0.034, 0.022) and / = (0.032, 0.032, 0.021), 
respectively. 

In the optically-thin approximation case, the effective 
cross-section is, by definition, zero at all times, and thus 
only the second method for calculation of ^ is applicable. 
By contrast with our results which took proper account of 
the optical depth in bound- free opacity, for which f <^ 1, our 
optically-thin results for ^ by the second approach yield / ~ 
1 [consistent with t he result obtained in this optica lly-thin 
approximation by iHaiman. Abel fc Madaul i200 J) 1. This 
enormous overestimate of ^ by the opticaUy-thin approxi- 
mation indicates that the optically-thin approximation used 
by previous authors is completely inadequate for determin- 
ing C- 

The reason for this significant discrepancy is easily un- 
derstood if we consider where and when most of the recom- 
binations take place. In Figure 1^ we show the evolution of 
^ for our illustrative simulations. During the initial R-type 
phase, the I-front deceleration phase, a significant fraction 
of the mass (~ 40% in our sample simulations) becomes 
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Figure 30. Evolution of tiie consumption of ionizing photons as 
given by equation l3Ui for all three spectra, as labelled. 



ionized almost instantly, just as does the entire minihalo 
in the optically-thin approximation. This gas is in the halo 
ourskirts, however, where the density is low and that den- 
sity drops even further as the gas expands in the evaporative 
outflow, into the IGM, making recombination in this gas in- 
efficient within 10 — 20 Myr. Most of the minihalo gas is 
initially shielded, however, and remains neutral for a much 
longer time. From this point on, the additional ionization 
of gas occurs only where the photoevaporation process is 
removing it by expelling it in a supersonic wind. For this 
gas, the density becomes too low too soon for significant re- 
combinations to occur there. Therefore, during most of the 
evolution, the bulk of the recombinations occur in a thin 
dense layer in the immediate vicinity of the I-front. Each 
gas atom spends only a short time in this layer, however, as 
it rushes supersonically away from the neutral region into 
the low-density expanding wind, hence experiencing many 
fewer recombinations than the optically-thin approximation 
would predict. 

As Figure 1301 shows, ^ increases gradually throughout 
the evaporation time, and there are noticeable differences 
among the three cases. For the hard QSO spectrum, the 
transition layer is thicker and penetrates deeper into the 
denser and colder parts of the halo, thus increasing the num- 
ber of recombinations per atom per evaporation time. How- 
ever, this same pre-heating effect shortens the evaporation 
time in this case, ultimately leading to a rough cancellation 
of the two effects and the same total ^ as in the BB 5e4 
case. The BB le5 spectrum is not as hard as in the QSO 
case, which decreases the number of recombinations in the I- 
front layer. The evaporation time in this case is also shorter 
than in the BB 5e4 case, due to a higher photoionization 



temperature and thus to faster outflow speed. Additionally, 
the higher temperatures in the wind for the BB le5 case 
(~ 10* K, as opposed to ~ 5000 K in the other two cases, 
see Figs. I24l and 1251 also leads to lower rates of recombina- 
tion in the wind and again lower ^. Thus, overall, the Pop. 
Ill sources appear significantly more efficient than Pop. 11 
or QSO sources in terms of the total number of ionizing 
photons needed to complete the photoevaporation process. 

In principle, the optically-thin approximation might be- 
come correct in the limit of xs ^ Xs.crit (see ii l2.3ll . in which 
case the 1-front is not trapped, the minihalo is ionized al- 
most instantly compared to the evaporation time, and the 
radiative transfer effects are less significant. However, for 
most minihalos this is not a relevant limit, since it requires 
excessively high ionizing flux of Fo > 100 — 1000 (depending 
on the halo mass and redshift of collapse) (see S 12.211 . 



5.9 Observational diagnostics 

5.9.1 Absorption lines 

Some observational signatures of the photoevaporation pro- 
cess are shown in Figures miliMI In Figures |^ and EH we 
show histograms of the column densities of H 1, He 1 and 
II, and C IV for minihalo gas at different radial velocities as 
seen along the symmetry axis at different times, to illustrate 
the kind of absorption lines which a photoevaporating mini- 
halo might produce. At early times, the minihalo gas resem- 
bles a weak damped Lya ("DLA") absorber with small ve- 
locity width (>10kms"^) and AfHi>10^°cm"^, with a Lya- 
Forest ("LF")-like red wing (velocity width > lOkms"^) 
with A^'hi ~ lO^^cm"^ on the side moving toward the source. 
As photoevaporation proceeds, this red wing increases in H 1 
column density to A'hi ~ lO^^cm"^. The He I profile mim- 
ics that of H 1 but with Nn^i/Nni ~ [He]/[H], and there 
is a weak C IV feature with New ~ 10^° (lO^^jr/cm-^ for 
the BB 5e4 (QSO, BB le5) cases, respectively, displaced 
in this same asymmetric way to the red side of the veloc- 
ity of peak H I column density, where r) = [C]/[C]0 x 10'^. 
For He II at early times, the QSO and BB 5e4 cases have 
A^'hcII ~ lO^^cm"^ at velocities close to those of the H 1 
peak, and a red wing shifted by ^ 10 km/sec to the red, with 
A^Hoii ~ lO'^^cm"^ which increases over time to 10^*cm~^. 
He 11 qualitatively follows the H I and He I profiles in these 
cases. For the BB le5 case, however, A'neii ~ lO^^cm"'^ in 
the HI peak, at early times, decreasing to iVHcii ~ lO^^cm"^ 
over time, with a red wing with A^hoII ~ lO^^cm"^ which 
increases to A^'hoII ~ 10^*cm~^. In this case, the He II profile 
notably diverges from the H I profile. After about an evapo- 
ration time, both H I and He I column densities decrease to 
LF-like values of lO" - 10^* cm , and shortly thereafter the 
minihalo is completely evaporated. Finally, the C IV column 
density is greatest in the red wing formed by the evaporative 
wind, except at the earliest times. 
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Figure 31. Observational diagnostics I. Absorption lines. Minihalo column densities (cm~ ) along symmetry axis for gas at different 
velocities, for photoevaporating minihalo. (a) (left) BB 5e4 case; (b) (right) QSO case. (Top) H I; (Middle) He I (solid) and He H 
(dotted); (Bottom) C IV (if [C]/[C]q = T] x 10~'^, then plotted values are Nqiv/^)- Each box labeled with time (in Myr) since arrival 
of intergalactic I-front. 



5.9.2 Ionization structure of metals 

Figures 1^ and 1^ show the spatial variation of the relative 
abundances of C, N, and O ions along the symmetry axis 
at t = 60 Myrs. While the QSO case shows the presence 
at 60 Myrs of low as well as high ionization stages of the 
metals, up to C V, N VI and O VI on the source side of 
the minihalo and GUI on the neutral side, the softer spec- 
trum of the BB 5e4 case yields less highly ionized gas both 
on the ionized side of the I-front (C III, N III, O III) and 
the neutral side (C II, N I, O I and II). The BB le5 case 
is intermediate between the other two cases, ionizing car- 
bon up to C V, but ionizing nitrogen and oxygen mostly to 
N IV and O IV, respectively, with smaller fractions of N V 
and O V and much smaller, although not negligible fractions 
of N VI and O VI than in the QSO case. In the optically- 
thin approximation case in Figure 1^ (right panels), similar 
high ionization stages are reached for the ionized gas ex- 
posed to the same photoionizing spectrum as in the results 
above which took proper account of optical depth, but the 
low ionization stages of neutral or partially ionized gas are 
entirely missing. In the optically-thin results, all the gas is 
ionized instantly, no neutral region of minihalo or shadow 
exists and the flow expands symmetrically in all directions, 
which would lead to quite different spectral signatures as 
compared to the realistic optically-thick cases shown above. 



6 SUMMARY AND CONCLUSIONS 

We have presented the first numerical gas dynamics and 
radiative transfer simulations of the encounter between an 
intergalactic I-front and a minihalo during cosmic reioniza- 
tion at high redshift, resulting in the photoevaporation of 
the minihalo gas. We have studied all stages of the photoe- 
vaporation process in detail, starting from the propagation 



of the fast, weak, R-type I-front in the low-density IGM, 
its structure and extent in pressure, temperature and ion- 
ization of H, He and a possible trace of metals for different 
possible source spectra, corresponding to quasars. Pop. Ill 
and Pop. II stars. We have demonstrated for the first time 
the phenomenon of I-front trapping when the global R-type 
I-front runs into a minihalo along its path and slows down 
inside the minihalo to an R-critical front, after which the 
front inside the minihalo transform from an R-type to a D- 
type front. We have shown that the heated and ionized gas 
behind the front then expands supersonically into the IGM, 
sweeping it up and shock-heating it. This process exposes 
further layers of the minihalo gas to the ionizing radiation 
from the external source, heating and evaporating them un- 
til all the gas is flnally expelled from the minihalo. For our 
illustrative case of a lO^M© halo, overtaken at z = 9 by 
an intergalactic I-front created by a source of 10^® ionizing 
photons per second at distance of 1 Mpc (or equivalently, 
a source of lO^'^sec"^ at 10 kpc), this process is completed 
within ~ 100 — 150 Myr, depending on the source spectrum. 

The results presented here are fully consistent with our 
earlier simulations of minihalo photoevaporation, going back 
to Shapiro, Raga, & Mellema (1997, 1998) and continu- 
ing through those in Shapiro & Raga (2000a, b; 2001), and 
Shapiro (2001), although there are some quantitative differ- 
ences which result from improvements we have made since 
then. However, these results differ from those which have 
appeared in the meantime in the rece nt literatu re. As men- 
tioned in §1.3, for example, .Barkana fc Loebl 111999.) con- 
sidered the photoevaporation of minihalos without gas dy- 
namics and without treating the propagation of the I-front. 
Their calculations are similar in spirit to (albeit more so- 
phisticated than) our analytical ISL approximation in H2.3I 
The result of their static approximation is that the larger 
minihalos (^ few x 10^ — 10® M©, depending on the red- 
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Figure 32. Same as Figure rm but for BB le5 case. 



shift and the source spectrum) were predicted not to evap- 
orate completely but rather to leave a significant fraction 
of the gas (20%-40% for QSO-type power-law spectrum for 
lO^M© halo at 2; = 8 — 20) still gravitationally bound to 
the minihalo. Our illustrative simulations demonstrate that 
when dynamical evolution is included this is not the correct 
outcome and all the gas is evaporated from the minihalos, 
leaving behind only a dark matter halo almost completely 
devoid of gas. In a companion paper, we will extend these 
studies to show that this conclusion holds for all minihalos, 
independent of their mass and redshift of collapse. 

We have also shown here that a hydrodynamical treat- 
ment of the photoevaporation of minihalos which fails to 
include radiative transfer (i.e. the zero optical depth ap- 
proximation) is also not adequate for describing this prob- 
lem. Such an approximation yields a spherically-symmetric 
outflow and uniformly high ionization structure which differ 
greatly from the more realistic results reported here which 
take radiative transfer into account. As a consequence, any 
observable features, such as absorption lines, would not be 
predicted correctly by this approximation. The difference 
between more realistic simulations like ours, which include 
radiative transfer, and those which neglect optical depth can 
have profound consequences for the theory of cosmic reion- 
ization, as welL 

Haiman. Abel fc MadaiJ (|200Jl pointed out the poten- 
tial importance of minihalos as sinks of ionizing photons 
during reionization, due to the increased recombination rate 
inside these high-density regions. In §2.5, we described how 
these authors estimated, in the optically-thin limit, the num- 
ber of ionizing photons per atom ^ = n-^/na needed to evap- 
orate a minihalo. They concluded from this that the clump- 
ing due to minihalos can easily raise the number of photons 
required to complete reionization by an order of magnitude 
or more compared to estimates that ignored the minihalos. 
This prompted them to identify a photon budget problem 
when comparing the requirements for reionizing the universe 



with the available photon supply. In this paper, we have used 
our self-consistent gas dynamics and radiative transfer sim- 
ulations to calculate ^ under more realistic assumptions. We 
have shown that the effect of ignoring radiative transfer and 
its feedback on the gas dynamics is to significantly overesti- 
mate ^. For the illustrative cases shown here, we find that ^ 
was overestimated by a factor as large as 30 — 50, depending 
on the ionizing spectrum. In terms of the efficiency factor 
in equations 1311 and I32I I. this means that we have found 
values of / ~ 1/30 to 1/50, depending on the spectrum of 
the ionizing sources. Optically-thin estimates also do not ac- 
count for differences due to differences in the external source 
spectra. We have shown that photoevaporation by Pop. Ill 
stars is significantly more efficient in terms of the net pho- 
ton consumption by this process than is photoevaporation 
by QSO or Pop. II stellar sources, even when all sources 
emit the same number of H ionizing photons per unit time. 
These eff ects may help alleviate the phot on budget problem 
posed bv lHaiman. Abel fc Madaul ll200j) . especially ff Pop. 
Ill sources dominated reionization. 

Although we find here that ^ is smaller than previously 
estimated, the effects of minihalos as screens and of their 
evaporation as a sink of ionizing photons during reionization 
may still be considerable. For one, we find ^ ~ 5, which is 
still large compared to unity. In addition, this value applies 
only to the particular illustrative case of a IO^Mq minihalo 
exposed to a source with fiux _Fo = 1 at z = 9. Only after we 
also calculate the values of ^ per minihalo for the full range of 
minihalo masses and redshifts and of source fluxes which can 
occur during reionization can we determine the full impact of 
the minihalos on that process. We will report those results 
in a companion paper to follow. However, we can already 
see that, if the reionization epoch was extended in time, 
as described in §1.2, with the final overlap of isolated H II 
regions occurring at 2 < 9, then the neutral regions which 
were ionized during this late phase had a significant fraction 
(> 30%) of their baryons collapsed into minihalos. In that 
case, minihalo photoevaporation by the global I-fronts which 
reionized the universe is likely to have dominated this final 
"overlap" phase. 

Apart from the destructive effects of intergalactic I- 
fronts on minihalo gas content described here, it has also 
been suggested that the same I-fronts might, instead, have 
a positive feedback effect on the minihalos rather than pho- 
toevaporating them entirely. It has been proposed, for ex- 
ample, that under certain conditions, external ionizing ra- 
diation can cause the implosion of the minihalo gas, leading 
to the formation of globular clusters ( Cen 200 L). We have 
not observed such an effect in our simulations. 

Finally, we note that intervening minihalos are ex- 
pected to be ubiquitous along the line of sight to high red- 
shift sources (Shapiro 2001). With photoevaporation times 
iov ^ 100 Myr, this process can continue down to redshifts 
significantly below z = 10. For stellar sources in the ACDM 
model, our simulations show that photoevaporation of a 
IO^Mq minihalo which begins at Zinitiai = 9 can take 150 
Myr to finish, at Zfinai = 7, during which time such minihalos 
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Figure 33. Observational diagnostics II: ionization structure of metals. C, N, and O ionic fractions along symmetry axis at f = 60 Myr, 
for photoevaporating minihalo (a) (left) BB 5e4 case; (b) (right) QSO case. 
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Figure 34. (a) (left) Same as Figure EH] but for BB le5 case, (b) (right) Same as Figure l33l but for optically-thin BB 5e4 case. 



can survive without merging into larger halos. Observations 
of the absorption spectra of high redshift sources hke those 
which reionized the universe should reveal the presence of 
these photoevaporative flows and provide a useful diagnostic 
of the reionization process. 
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APPENDIX A: A TIME-DEPENDENT MODEL 
FOR VIRIALIZED HALOS WITH 
COSMOLOGICAL INFALL 

The initial conditions for our simulations are given by the 
TIS profile for both the gas and the dark matter density 
for the virialized halo, wit h zero bulk velocity, as described 
by llliev fc Shapirol i200 J) , and a matching spherical, self- 
similar infall profile as described bv lBertschingeJ il985l) for 
the density and infall velocity outside the halo, appropri- 
ately generalized to the case of a low-density background 
universe models (Iliev & Shapiro 2001, Appendix A). The 
match of the post-shock gas in this infall solution to the 
TIS is discussed in detail in § 7.2 of IShapiro. Iliev fc Ragal 
1^9^. This match is made possible by the fact that the TIS 
radius rt is almost coincident with the shock radius of the 
self-similar infall solution rs, with the latter larger than the 
former by only 1.8%.* In order to match the two solutions 
seamlessly, we continue both the density and the velocity 
profiles of the infall solution down to rt- For the velocity 
profile outside the TIS outer radius we use the pre-shock 
infall pro file of the sel f -simil ar solution. Following the no- 
tation of iBertschingeil il985r) . we define the dimensionless 
radius A = r/rta(i). Here rta is the turnaround radius of the 
shell which is just reaching maximum expansion at time t, 
whose time-dependence is given by rta(t) = rt!i(ti)(t/ti)^^^ , 
where ti is our initial time. Due to self-similarity of the so- 
lution every feature in radius occurs at some fixed value 
of A. For example, the location of the shock is given by 
As = 0.338976, and the time- varying shock radius must fol- 
low rs{t) = rta(t)As. Similarly, the radial parameters of the 
TIS profile, like core radius ro and Vt, follow the same time 
dependence with A = (Ao,At) = (0.0113175,0.3327339), re- 
spectively. 

The physical velocity i; at a given radius r ^ rt t a. time 
t is given by the following parametric solution: 



* In this Appendix, spherical symmetry applies throughout, so 
we use the variable "r" to mean the spherical radius here. 



X{e) = sm - 
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and 



t (l-cos( 



yfl ^ ro(t) 
|2 Xot 



(Al) 



(A2) 



where ^ ^ ^ 2n. For the TIS density profile we use the 
fitting for mula to the exact nume rical results given in Ap- 
pendix B of llliev fc Shapirol i200lf) . The resulting initial con- 
dition for our simulations is shown in Figure 2] 

In order to evolve the initial dark matter profile with 
time self-consistently, we assume that the mass of the TIS 
grows over time as a result of the infall, and our initial so- 
lution for the TIS gives way to a sequence of equilibrium 
TIS models, where the TIS is specified at each time by the 
mass M{t) and redshift Zcoii at that time. Below we give the 
time- dependences of all relevant quantities. 

The time-dependence of M{t) corresponds to that of the 
sphere bounded by the accretion shock in the Bertschinger 
solution. The central density po of our TIS solution is just 
proportional to the mean background density at the same 
epoch. At the early times of interest to us here, this depen- 
dence is simply po{t) — po{ti){t/ti)~'^ . With this scaling we 
can write the growth of the total mass with time as 



MTis(i) = MTis(ii 



Mtis 



rUt) 



Pb{t) 



pbiti 



t\^/'i 



(A3) 



In fact, for any radius which grows in time self-similarly (i.e. 
keeping its A- value constant) the mass interior to this radius 
must grow according to 



m{X,t) = m{X,ti) J 



t \2/3 



(A4) 



The gravitational acceleration, g, at radius r of this as- 
sumed matter distribution is given simply in terms of the 
mass m(r) interior to this radius: g = Gm{r)/r'^ . For fixed 
A, m{X,t) (X t'^^^, so we need only determine m{r,t) = 
m[X{r,t),t] at fixed r for different times. This is given by 



m{X,t) = M(A) 
where 



47rp6(ii)r| fty/^ 



3A| 



X{r,t) = 



Ao 



ro{ti) 



») iti) 



-8/9 



(A5) 



(A6) 



(A) For A ^ At = 0.332734, the infall mass profile is 



M(A) 



9 



A^D(A) 



V{X) - |a 



where V^(A) was defined in equation l|A2|l above and 

^ 9ie-sm9f 

^ ' 16sin«(6»/2){l + 3[l-3V(A)/(2A)]} 

jBertschingeilll985h . 



(A7) 



(A8) 
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(B) For A ^ At the TIS mass profile is derived by enforc- 
ing liydrostatic equilibrium inside the TIS, using the density 
profile fitting formula and the isothermal virial temperature 
to solve for the gravitational acceleration needed to balance 
the pressure force: 



where (^,a^B,fe^) = (21.38,9.08,19.81,14.62), ^ = r/ro, 
ro is the TIS core radius, and pb is the mean background 
density. 




X 



[A/(a2 5/(62 



(A9) 
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